trajectory.qmd

---
title: "Trajectories"
---

### Standard analysis

We load our usual example data

```{r}
suppressPackageStartupMessages({
  library( tidyverse )
  library( Matrix )
  library( sparseMatrixStats )
  library( Seurat ) })

ReadMtx( "~/Downloads/ifnagrko/ifnagrko_raw_counts.mtx.gz",
    "~/Downloads/ifnagrko/ifnagrko_obs.csv",
    "~/Downloads/ifnagrko/ifnagrko_var.csv",
    cell.sep=",", feature.sep=",", skip.cell=1, skip.feature=1, 
    mtx.transpose=TRUE) -> count_matrix
```

```{r}
count_matrix %>%
CreateSeuratObject() %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA( npcs=20 ) %>%
FindNeighbors( dims=1:20 ) %>%
FindClusters( resolution=0.5 ) %>%
RunUMAP( dims=1:20 ) -> seu
```

```{r}
UMAPPlot( seu, label=TRUE ) + coord_equal()
```

This time, we will concetrate on the long elongated main structure, which we
will call the "lineage" in the following. It is a snapshot of the development
of astrocytes that act as neural stem cells and become transient amplifying 
progenitors (TAPs) which undergo cell cycle, i.e., divide and multiply, and the 
turn into neuroblasts and finally neurons.

To orient us in the plot, we highlight the expression of Aqp4 (aquaporin-4, a marker for 
astrocytes), Mki67 (a marker for prliferating, i.e., dividing cells), Dcx 
(doublecortin, a marker for neuroblasts) and Gria1 (Glutamate ionotropic receptor, 
AMPA type, subunit 1; a marker for mature neurons)

```{r}
FeaturePlot( seu, c( "Aqp4", "Mki67", "Dcx", "Gria1" ) )
```

We conclude that the lineage is well covered by the following clusters

```{r}
lineage_clusters <-  c( 10, 9, 0, 13, 14, 3, 5, 6, 11, 1, 2, 7 ) 
```

Just out of curiosity, we also try to identify clusters 4 and 8:

```{r}
presto::wilcoxauc(
  LayerData(seu),
  factor(case_when(
    seu$seurat_clusters %in% lineage_clusters ~ "lineage",
    seu$seurat_clusters %in% c( 4, 8 ) ~ "4_or_8",
    TRUE ~ "other" )) ) -> wa

head(wa)
```

```{r}
presto::top_markers( wa )
```

As a try, I've asked ChatGPT what these genes point to and it [replied](https://chatgpt.com/share/f2a5e86b-9909-4c50-8b6b-798f9596df71) that the cells in our clusters 4 and 8 are oligondendrocyte precursor cells (OPCs). This matches
my expectation.

### Aim: Trajectory

Our first aim for today is to fit a "pseudotime trejectory" to the lineage. This means
thatw e want to assign to each cell in the lineage a real number, which we call
it "pseudotime" that monotonally increases along the putative developmental trajectory 
from astrocytic neural stem cells via TAPs and neuroblasts to neurons.

We will do this by fitting a "principal curve", i.e. a curve in PCA space that tracks
along the lineage and is fitted such that the squared sum of the cells' distance to
their respectively closest point on the curve (their projection image) is minimal.
The distance of this projection image to the curve start (measured along the curve)
will be used as pseudotime.

The "principal curve" method is described in detail in [this section of the lecture notes](principal_curves.html).

Here, we will use the function from the `princurve` package.

As preparation, we first explore distances with sleepwalk:

```{r eval=FALSE}
sleepwalk::sleepwalk( Embeddings(seu,"umap"), Embeddings(seu,"pca") )
```
We notice that the cycling cells have a lot of distance to the non-cycling
lineage cells. This will cause problems because the principal curve alorithm
cannot deal with loops, i.e., the curve should pass besides the cell-cycle loop.
However, the distances to the curve will then become large right in the middle, deflecting
the curve.

Therefore, let's exclude clusters 6 and 11:

```{r}
lineage_clusters_2 <- setdiff( lineage_clusters, c( 6, 11 ) )
```

### Fitting the principal curve

Now, we could the `princomp` package to fit the principal curve

```{r}
princurve::principal_curve( 
  Embeddings(seu,"pca")[ seu$seurat_clusters %in% lineage_clusters_2, ], 
  df=10, trace=TRUE, approx_points=1000 ) -> prc
```

This function returns for each cell a pseudotime value `lambda` and a projected position on
the curve in PCA space, `s`.

```{r}
Embeddings(seu,"umap") %>%
as_tibble( rownames="cell" ) %>%
left_join( enframe( prc$lambda, "cell", "lambda" ) ) %>%
ggplot +
  geom_point( aes( x=umap_1, y=umap_2, col=lambda ), size=.3 ) +
  coord_equal() + scale_color_viridis_c(option="D")
```

We can assign a pseudotime to the remaining cells by finding the closest curve point:

```{r}
FNN::get.knnx( prc$s, Embeddings(seu,"pca"), 1 ) -> nnres

prc$lambda[ nnres$nn.index[,1] ] %>%
  { ( max(.) - . ) / max(.) } %>%
  set_names( rownames(Embeddings(seu,"pca")) ) -> seu$pt
```

Her we have rescaled the pseudotime to [0;1] and also reversed the dirction,
so thjat it now increases from stem cells towards neuroblasts.

```{r}
Embeddings(seu,"umap") %>%
as_tibble( rownames="cell" ) %>%
left_join( enframe( seu$pt, "cell", "pt" ) ) %>%
ggplot +
  geom_point( aes( x=umap_1, y=umap_2, col=pt ), size=.3 ) +
  coord_equal() + scale_color_viridis_c(option="D")
```

We should also check how for each cell is from the curve

```{r}
nnres$nn.dist[,1] %>%
  set_names( rownames(Embeddings(seu,"pca")) ) -> seu$dist_to_curve
  
Embeddings(seu,"umap") %>%
as_tibble( rownames="cell" ) %>%
left_join( enframe( seu$dist_to_curve, "cell", "dist" ) ) %>%
ggplot +
  geom_point( aes( x=umap_1, y=umap_2, col=dist ), size=.3 ) +
  coord_equal() + scale_color_viridis_c( option="D", trans="log10", direction=-1 )  
```

### Expression dynamics

Here is a plot showing the expression of one gene, Slc1a3 (Glast), along the 
pseudotime:

```{r}
tibble(
  pt = seu$pt, 
  dist = seu$dist_to_curve,
  in_lineage = seu$seurat_clusters %in% lineage_clusters,
  expr = LayerData(seu)["Slc1a3",] ) %>%
mutate( expr = ifelse( expr>0, expr, runif( n(), -.2, 0 ) ) ) -> tbl

tbl %>%
filter( in_lineage ) %>%
ggplot +
  geom_point( aes( x=pt, y=expr, col=dist ), size=.3 ) +
  scale_color_viridis_c()
```
In order to make more apparent how many point are on the zero line, this line has been broadened.

We now fit a smooth curve through this scatter plot. As we want to do this properly,
we use locfit with Poisson GLM, i.e., wo don't use the log-normalized
values but the raw counts. (Details [here](smoothing.html).)

First we assemble the data

```{r}
library( locfit )

tibble(
  pt = seu$pt, 
  dist = seu$dist_to_curve,
  in_lineage = seu$seurat_clusters %in% lineage_clusters,
  count = LayerData(seu,"count")["Slc1a3",],
  total = colSums( LayerData(seu,"count") ),
  expr = LayerData(seu)["Slc1a3",] ) %>%
mutate( expr = ifelse( expr>0, expr, runif( n(), -.2, 0 ) ) ) -> tbl

head(tbl)
```

The we run `lucfit`:
```{r}
fit <- locfit( count ~ pt, tbl, weight=total, family="poisson" )

fit
```

We evaluate the fitted curve along a value grid:

```{r}
tibble( pt = seq( 0, 1, length.out=1000 ) ) %>%
mutate( y = predict( fit, pt ) ) -> tbl_fit

head(tbl_fit)
```

Now we can do the plot:

```{r}
tbl %>%
filter( in_lineage ) %>%
ggplot( aes( x=pt ) ) +
  geom_point( aes( y = count/total + 1e-4, col=dist ), size=.3 ) +
  geom_line( aes( y = y ), data=tbl_fit, col="magenta" ) +
  scale_color_viridis_c() + scale_y_log10()
```

#### Many genes

We can run this for several genes. We pick the 10 genes with the highest variance
of expression along the lineage:

```{r}
LayerData(seu)[ , seu$seurat_clusters %in% lineage_clusters ] %>%
rowVars() %>%
sort( decreasing=TRUE ) %>%
head(10) %>% names() -> genes

tg2 <- seq( 0, 1, length.out=300 )

sapply( genes, function(gene) {
  fit <- locfit( LayerData(seu,"count")[gene,] ~ seu$pt, 
          weight=seu$nCount_RNA, family="poisson" )
  cat(".") 
  predict( fit, tg2 )
} ) %>% t() -> fits

fits[1:5,1:5]
```

This time, we have evaluated the smoothed curve at a grid of 300 values and got a 
matrix with one row per gene and one column for each of the 300 time points.

To visualize this, we use a heatmap:

```{r}
image( t(fits) )
```

The genes have different dynamic range. Hence, we should divide each row by its maximum:

```{r}
fitsz <- fits / rowMaxs(fits)
image( t(fitsz), yaxt="n" )
axis( 2, seq( 0, 1, length.out=nrow(fitsz) ), rownames(fitsz), las=2, cex.axis=.5 )
```

Now, let's also sort the rows by the positions of these maxima, and replace the colour scale:

```{r}
fitszs <- fitsz[ order( -apply( fitsz, 1, which.max ) ), ]
image( t(fitszs), yaxt="n", col=viridisLite::viridis(300) )
axis( 2, seq( 0, 1, length.out=nrow(fitszs) ), rownames(fitszs), las=2, cex.axis=.5 )
```

The same now with the top hundred genes:

```{r}
LayerData(seu)[ , seu$seurat_clusters %in% lineage_clusters ] %>%
rowVars() %>%
sort( decreasing=TRUE ) %>%
head(100) %>% names() -> genes

sapply( genes, function(gene) {
  fit <- locfit( LayerData(seu,"count")[gene,] ~ seu$pt, 
          weight=seu$nCount_RNA, family="poisson" )
  cat(".") 
  predict( fit, tg2 )
} ) %>% t() -> fits
```

```{r fig.width=6,fig.height=15}
fitsz <- fits / rowMaxs(fits)
fitszs <- fitsz[ order( -apply( fitsz, 1, which.max ) ), ]
image( t(fitszs), yaxt="n", col=viridisLite::viridis(300) )
axis( 2, seq( 0, 1, length.out=nrow(fitszs) ), rownames(fitszs), las=2, cex.axis=.5 )
```