Chromaffin differentiation analysis

The notebook shows an example of basic velocity anlaysis of E12.5 chromaffin data, which is a SMART-seq2 dataset. It shows how to estimate gene-relative velocity (with different pooling options), as well as how to estimate using alternative gene-structure based model.

Data pre-processing

We assume that the SMART-seq2 data has been demultiplexed and aligned, so that we have one bam file per cell. The recommended processing step is then to use velocyto.py command line tool to annotate spliced, unspliced and spanning reads in the measured cells:

velocyto run_smartseq2 -d 1 `find ./ data/e12.5.bams/ -name "*.bam" -print` mm10.genes.gtf

If you want to re-run this command, please download bam files (and genes.refFlat) from here, and extract it (tar xvf bams.tar) in the working directory. note: the file is fairly large - 5.6 GB! Running this command will produce a velocyto directory with two files:

1. onefilepercell_A1_unique_and_others_J2CH1.loom - a loom file containing spliced, unspliced and spanned matrices
2. dump/onefilepercell_A1_unique_and_others_J2CH1.hdf5 - an HDF5 file containing detailed molecular mapping information which we will use to estimate parameters for the gene-structure based model. You don’t need this for regular velocity estimation, and in this case -d 1 option may be skipped when running velocyto.py

note: if you don’t want to re-run the command, you can download the velocyto.py output folder here

Data loading

Load the velocyto package:

library(velocyto.R)

Load the data and genome annotations:

Read in count matrices from the loom file:

ldat <- read.loom.matrices("data/velocyto/onefilepercell_A1_unique_and_others_J2CH1.loom")

(instead, for the purposes of the example we’ll just read in the resulting structure from the rds file)

ldat <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/ldat.rds"))
str(ldat)

List of 4
 $ spliced  :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:2724270] 3 9 10 13 14 17 19 20 28 29 ...
  .. ..@ p       : int [1:385] 0 7100 14272 21376 28491 36017 43375 50956 57637 64998 ...
  .. ..@ Dim     : int [1:2] 23420 384
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:384] "onefilepercell_A1_unique_and_others_J2CH1:A1_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A10_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A11_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A12_unique.bam" ...
  .. ..@ x       : num [1:2724270] 24 3 22 136 48 21 15 32 109 13 ...
  .. ..@ factors : list()
 $ unspliced:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:1528882] 0 4 10 13 14 17 19 20 23 29 ...
  .. ..@ p       : int [1:385] 0 3911 8043 12057 16176 20795 25046 29389 33438 37382 ...
  .. ..@ Dim     : int [1:2] 23420 384
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:384] "onefilepercell_A1_unique_and_others_J2CH1:A1_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A10_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A11_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A12_unique.bam" ...
  .. ..@ x       : num [1:1528882] 5 23 4 6 8 2 70 14 1 15 ...
  .. ..@ factors : list()
 $ ambiguous:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:198206] 3 19 30 121 144 145 188 295 366 387 ...
  .. ..@ p       : int [1:385] 0 499 993 1562 2070 2639 3184 3785 4321 4821 ...
  .. ..@ Dim     : int [1:2] 23420 384
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:384] "onefilepercell_A1_unique_and_others_J2CH1:A1_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A10_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A11_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A12_unique.bam" ...
  .. ..@ x       : num [1:198206] 17 1 14 24 9 1 5 23 1 2 ...
  .. ..@ factors : list()
 $ spanning :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:443099] 3 17 40 80 111 132 146 161 200 216 ...
  .. ..@ p       : int [1:385] 0 1089 2250 3556 4768 6201 7494 8878 10237 11349 ...
  .. ..@ Dim     : int [1:2] 23420 384
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  .. .. ..$ : chr [1:384] "onefilepercell_A1_unique_and_others_J2CH1:A1_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A10_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A11_unique.bam" "onefilepercell_A1_unique_and_others_J2CH1:A12_unique.bam" ...
  .. ..@ x       : num [1:443099] 1 2 1 3 2 1 1 1 1 1 ...
  .. ..@ factors : list()

Let’s reduce the cell names to the short well labels:

ldat <- lapply(ldat,function(x) {
  colnames(x) <-  gsub("_unique.bam","",gsub(".*:","",colnames(x)))
  x
})

Read in cell cluster assignment and tSNE embedding used in the Furlan et al. (Science’17).

cell.colors <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/cell.colors.rds"))
emb <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/embedding.rds"))

Gene filtering

Spliced expression magnitude distribution across genes:

hist(log10(rowSums(ldat$spliced)+1),col='wheat',xlab='log10[ number of reads + 1]',main='number of reads per gene')

Set up expression matrices, filtering genes to leave those that exceed some pre-defined g to the average expression magnitude

# exonic read (spliced) expression matrix
emat <- ldat$spliced;
# intronic read (unspliced) expression matrix
nmat <- ldat$unspliced
# spanning read (intron+exon) expression matrix
smat <- ldat$spanning;
# filter expression matrices based on some minimum max-cluster averages
emat <- filter.genes.by.cluster.expression(emat,cell.colors,min.max.cluster.average = 5)
nmat <- filter.genes.by.cluster.expression(nmat,cell.colors,min.max.cluster.average = 1)
smat <- filter.genes.by.cluster.expression(smat,cell.colors,min.max.cluster.average = 0.5)
# look at the resulting gene set
length(intersect(rownames(emat),rownames(nmat)))

[1] 8548

# and if we use spanning reads (smat)
length(intersect(intersect(rownames(emat),rownames(nmat)),rownames(smat)))

[1] 1696

Several variants of velocity estimates using gene-relative model

We’ll start with what is perhaps the most robust estimate, that combines cell kNN pooling with the gamma fit based on an extreme quantiles:

Using min/max quantile fit, in which case gene-specific offsets do not require spanning read (smat) fit. Here the fit is based on the top/bottom 5% of cells (by spliced expression magnitude).

fit.quantile <- 0.05;
rvel.qf <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,kCells = 5,fit.quantile = fit.quantile)

calculating cell knn ... done
calculating convolved matrices ... done
fitting gamma coefficients ... done. succesfful fit for 8548 genes
filtered out 1306 out of 8548 genes due to low nmat-emat correlation
filtered out 754 out of 7242 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done

We visualize the velocities by projecting observed and extrapolated cells onto the first 5 PCs:

pca.velocity.plot(rvel.qf,nPcs=5,plot.cols=2,cell.colors=ac(cell.colors,alpha=0.7),cex=1.2,pcount=0.1,pc.multipliers=c(1,-1,-1,-1,-1))

log ... pca ... pc multipliers ... delta norm ... done
done

Fitting of individual genes can be visualized using “show.gene” option. To save time, we’ll pass previously-calculated velocity (rvel.qf) to save calculation time:

# define custom pallet for expression magnitude
gene.relative.velocity.estimates(emat,nmat, kCells = 5,fit.quantile = fit.quantile,old.fit=rvel.qf,show.gene='Chga',cell.emb=emb,cell.colors=cell.colors)

calculating convolved matrices ... done
[1] 1

Alternatively, wen ca use spanning reads (smat) to fit the gene offsets. This will result in more accurate offset estimates, but for much fewer genes (spanning reads are rare). Note that here we also use optional diagona.quantiles setting to estimate extreme quantiles on a normalized sum of splcied and unspliced signal:

rvel <- gene.relative.velocity.estimates(emat,nmat,smat=smat, kCells = 5, fit.quantile=fit.quantile, diagonal.quantiles = TRUE)

calculating cell knn ... done
calculating convolved matrices ... done
fitting smat-based offsets ... done
fitting gamma coefficients ... done. succesfful fit for 1696 genes
filtered out 26 out of 1696 genes due to low nmat-smat correlation
filtered out 138 out of 1670 genes due to low nmat-emat correlation
filtered out 14 out of 1532 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done

We can visualize the velocity in PCA space:

pca.velocity.plot(rvel,nPcs=5,plot.cols=2,cell.colors=ac(cell.colors,alpha=0.7),cex=1.2,pcount=0.1,pc.multipliers=c(1,-1,1,1,1))

log ... pca ... pc multipliers ... delta norm ... done
done

Here we calculate the most basic version of velocity estimates, using relative gamma fit, without cell kNN smoothing (i.e. actual single-cell velocity):

rvel1 <- gene.relative.velocity.estimates(emat,nmat,deltaT=1,deltaT2 = 1,kCells = 1, fit.quantile=fit.quantile)

fitting gamma coefficients ... done. succesfful fit for 8548 genes
filtered out 783 out of 8548 genes due to low nmat-emat correlation
filtered out 1330 out of 7765 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done

pca.velocity.plot(rvel1,nPcs=5,plot.cols=2,cell.colors=ac(cell.colors,alpha=0.7),cex=1.2,pcount=0.1,pc.multipliers=c(1,-1,1,1,1))

log ... pca ... pc multipliers ... delta norm ... done
done

Visualization on an existing embedding

Here we use t-SNE embedding from the original publication (in emb variable).

vel <- rvel; arrow.scale=3; cell.alpha=0.4; cell.cex=1; fig.height=4; fig.width=4.5;
show.velocity.on.embedding.cor(emb,vel,n=100,scale='sqrt',cell.colors=ac(cell.colors,alpha=cell.alpha),cex=cell.cex,arrow.scale=arrow.scale,arrow.lwd=1)

delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done

Alternatively, the same function can be used to calculate a velocity vector field:

show.velocity.on.embedding.cor(emb,vel,n=100,scale='sqrt',cell.colors=ac(cell.colors,alpha=cell.alpha),cex=cell.cex,arrow.scale=arrow.scale,show.grid.flow=TRUE,min.grid.cell.mass=0.5,grid.n=20,arrow.lwd=2)

delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done
grid estimates ... grid.sd= 1.731696  min.arrow.size= 0.03463392  max.grid.arrow.length= 0.09156871  done

Velocity estimate based on gene structure

To estimate velocity based on gene structure parameters, we need to parse the gtf file as well as the debug hdf5 output of velocyto.py (-d option) which contains per-exon mapping information. You can fetch the exact gtf file used here. First, however, we will compile information on the internal priming sites expected in this genome (mouse UCSC mm10 assembly):

require(BSgenome.Mmusculus.UCSC.mm10)
ip.mm10 <- find.ip.sites('data/genes.gtf',Mmusculus,'mm10')

This needs to be ran once per genome, and can be saved for the future. To skip calculations, we just load the pre-calculated internal priming info from an rds file:

ip.mm10 <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/ip.mm10.rds"))

Now we will read velocyto.py HDF5 output:

gene.info <- read.gene.mapping.info("data/velocyto/dump/onefilepercell_A1_unique_and_others_J2CH1.hdf5",internal.priming.info=ip.mm10,min.exon.count=10);

Will load pre-calculated result instead of evaluating above: This needs to be ran once per genome, and can be saved for the future. To skip calculations, we just load the pre-calculated internal priming info from an rds file:

gene.info <- readRDS(url("http://pklab.med.harvard.edu/velocyto/chromaffin/gene.info.rds"))

Genome-wide model fit:

# start with unfiltered matrices, as we can use more genes in these types of estimates
emat <- ldat$spliced; nmat <- ldat$unspliced; smat <- ldat$spanning
emat <- filter.genes.by.cluster.expression(emat,cell.colors,min.max.cluster.average = 7)
gvel <- global.velcoity.estimates(emat, nmat, rvel, base.df=gene.info$gene.df, smat=smat, deltaT=1, kCells=5, kGenes = 15, kGenes.trim = 5, min.gene.cells = 0, min.gene.conuts = 500)

filtered out 4 out of 8990 genes due to low emat levels
filtered out 1072 out of 8986 genes due to insufficient exonic or intronic lengths
filtered out 204 out of 7914 genes due to excessive nascent counts
using relative slopes for 1221 genes to fit structure-based model ... with internal priming info ... 80.5% deviance explained.
predicting gamma ... done
refitting offsets ... calculating cell knn ... done
calculating convolved matrices ... done
fitting smat-based offsets ... done
fitting gamma coefficients ... done. succesfful fit for 7694 genes
filtered out 1337 out of 7551 genes due to low nmat-smat correlation
filtered out 899 out of 6214 genes due to low nmat-emat correlation
filtered out 440 out of 5315 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done
re-estimated offsets for 6214 out of 7710 genes
calculating convolved matrices ... done
calculating gene knn ... done
estimating M values ... adjusting mval offsets ... re-estimating gamma ... done
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done

pca.velocity.plot(gvel,nPcs=5,plot.cols=2,cell.colors=ac(cell.colors,alpha=0.7),cex=1.2,pcount=0.1,pc.multipliers=c(1,-1,-1,1,1))

log ... pca ... pc multipliers ... delta norm ... done
done

We can visualize the two estimates side-by side on a pre-calculated (published) tSNE embedding:

par(mfrow=c(1,2), mar = c(2.5,2.5,2.5,1.5), mgp = c(2,0.65,0), cex = 0.85);
arrow.scale=3; cell.alpha=0.4; cell.cex=1; fig.height=4; fig.width=4.5;
#pdf(file='tsne.rvel_gvel.plots.pdf',height=6,width=12)
show.velocity.on.embedding.cor(emb,rvel,n=100,scale='sqrt',cell.colors=ac(cell.colors,alpha=cell.alpha),cex=cell.cex,arrow.scale=arrow.scale,arrow.lwd=1,main='gene-relative esitmate',do.par=F)

delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done

show.velocity.on.embedding.cor(emb,gvel,n=100,scale='sqrt',cell.colors=ac(cell.colors,alpha=cell.alpha),cex=cell.cex,arrow.scale=arrow.scale,arrow.lwd=1,main='gene-structure estimate',do.par=F)

delta projections ... sqrt knn ... transition probs ... done
calculating arrows ... done

#dev.off()

Or by jointly embedding both observed and extrapolated cells into the same tSNE space

#pdf(file='tsne.shift.plots.pdf',height=6,width=12)
par(mfrow=c(1,2), mar = c(2.5,2.5,2.5,1.5), mgp = c(2,0.65,0), cex = 0.85);
x <- tSNE.velocity.plot(rvel,nPcs=15,cell.colors=cell.colors,cex=0.9,perplexity=200,norm.nPcs=NA,pcount=0.1,scale='log',do.par=F,main='gene-relative estimate')

rescaling ... log ... pca ... delta norm ... tSNE ...Read the 768 x 15 data matrix successfully!
OpenMP is working...
Using no_dims = 2, perplexity = 200.000000, and theta = 0.500000
Computing input similarities...
Normalizing input...
Building tree...
Done in 1.46 seconds (sparsity = 0.903880)!
Learning embedding...
Iteration 50: error is 45.227904 (50 iterations in 1.37 seconds)
Iteration 100: error is 45.218709 (50 iterations in 1.42 seconds)
Iteration 150: error is 44.385833 (50 iterations in 1.29 seconds)
Iteration 200: error is 44.370694 (50 iterations in 1.60 seconds)
Iteration 250: error is 44.369845 (50 iterations in 1.70 seconds)
Iteration 300: error is 0.164786 (50 iterations in 1.22 seconds)
Iteration 350: error is 0.154593 (50 iterations in 1.20 seconds)
Iteration 400: error is 0.154770 (50 iterations in 1.17 seconds)
Iteration 450: error is 0.154346 (50 iterations in 1.17 seconds)
Iteration 500: error is 0.154502 (50 iterations in 1.18 seconds)
Iteration 550: error is 0.154393 (50 iterations in 1.19 seconds)
Iteration 600: error is 0.154297 (50 iterations in 1.18 seconds)
Iteration 650: error is 0.154556 (50 iterations in 1.16 seconds)
Iteration 700: error is 0.155218 (50 iterations in 1.17 seconds)
Iteration 750: error is 0.154701 (50 iterations in 1.18 seconds)
Iteration 800: error is 0.154848 (50 iterations in 1.20 seconds)
Iteration 850: error is 0.155505 (50 iterations in 1.21 seconds)
Iteration 900: error is 0.155638 (50 iterations in 1.20 seconds)
Iteration 950: error is 0.155412 (50 iterations in 1.21 seconds)
Iteration 1000: error is 0.155655 (50 iterations in 1.19 seconds)
Fitting performed in 25.22 seconds.
delta norm ... done

x <- tSNE.velocity.plot(gvel,nPcs=15,cell.colors=cell.colors,cex=0.9,perplexity=200,norm.nPcs=NA,pcount=0.1,scale='log',do.par=F,main='gene-structure estimate')

rescaling ... log ... pca ... delta norm ... tSNE ...Read the 768 x 15 data matrix successfully!
OpenMP is working...
Using no_dims = 2, perplexity = 200.000000, and theta = 0.500000
Computing input similarities...
Normalizing input...
Building tree...
Done in 1.29 seconds (sparsity = 0.905877)!
Learning embedding...
Iteration 50: error is 45.043481 (50 iterations in 1.33 seconds)
Iteration 100: error is 45.043481 (50 iterations in 1.32 seconds)
Iteration 150: error is 45.043481 (50 iterations in 1.25 seconds)
Iteration 200: error is 45.043481 (50 iterations in 1.27 seconds)
Iteration 250: error is 45.043481 (50 iterations in 1.26 seconds)
Iteration 300: error is 0.172803 (50 iterations in 1.17 seconds)
Iteration 350: error is 0.156917 (50 iterations in 1.28 seconds)
Iteration 400: error is 0.157191 (50 iterations in 1.50 seconds)
Iteration 450: error is 0.157446 (50 iterations in 1.17 seconds)
Iteration 500: error is 0.156811 (50 iterations in 1.11 seconds)
Iteration 550: error is 0.156514 (50 iterations in 1.27 seconds)
Iteration 600: error is 0.156040 (50 iterations in 1.12 seconds)
Iteration 650: error is 0.156383 (50 iterations in 1.11 seconds)
Iteration 700: error is 0.155723 (50 iterations in 1.08 seconds)
Iteration 750: error is 0.155851 (50 iterations in 0.89 seconds)
Iteration 800: error is 0.156389 (50 iterations in 1.06 seconds)
Iteration 850: error is 0.156084 (50 iterations in 1.04 seconds)
Iteration 900: error is 0.156497 (50 iterations in 1.06 seconds)
Iteration 950: error is 0.156577 (50 iterations in 1.01 seconds)
Iteration 1000: error is 0.157346 (50 iterations in 1.06 seconds)
Fitting performed in 23.33 seconds.
delta norm ... done

#dev.off()

Cell trajectory modeling

A similar function can be used to model central trajectories by directed diffusion on embedding. The main parameters are set up by sigma (which limits the range of how far a cell can jump in terms of distance) and n (how many nearest neighbors are being considered for jumps). The results are sensitive to these parameters, as we don’t have a good way of assessing how much the directional velocity component should compare with random Brownian motion of a cell with the manifold. For instance, relaxing (increasing) sigma, in particular will eventually lead to sympathoblast cells “jumping” the gap into the into the chromaffin differentiation part.

Warning: this simulation takes some time (e.g. a couple of minutes on 40 cores).

x <- show.velocity.on.embedding.eu(emb,rvel,n=40,scale='sqrt',cell.colors=ac(cell.colors,alpha=cell.alpha),cex=cell.cex,nPcs=30,sigma=2.5,show.trajectories=TRUE,diffusion.steps=400,n.trajectory.clusters=15,ntop.trajectories=1,embedding.knn=T,control.for.neighborhood.density=TRUE,n.cores=40) 

sqrt scale ... reducing to 30 PCs ... distance ... sigma= 2.5  beta= 1  transition probs ... embedding kNN ... done
simulating diffusion ... constructing path graph ... tracing shortest trajectories ... clustering ... done.