Chromaffin differentiation analysis

The notebook shows anlaysis of a SMART-seq2 dataset, where we start with per-cell bam files. It shows how to estimate gene-relative velocity (with different pooling options), as well as how to estimate gene-relative velocity.

Data loading

Load the velocyto package:

library(velocyto.R)

Load the data and genome annotations:

(this block will note evalaute by default, as it takes a while)

# if you want to run this, please download bam files (and genes.refFlat) from
# http://pklab.med.harvard.edu/velocyto/chromaffin/bams.tar
# and extract it ("tar xvf bams.tar") in the working directory.
# note: the file is fairly large - 5.6 GB! 
path <- "data/e12.5.bams"
files <- system(paste('find',path,'-name "*unique.bam" -print'),intern=T)
names(files) <- gsub(".*\\/(.*)_unique.bam","\\1",files)
# parse gene annotation, annotate bam file reads
dat <- read.smartseq2.bams(files,"data/genes.refFlat",n.cores=40)

(instead, we read in the resulting structure from the rds file)

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

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(dat$emat)+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 <- dat$emat;
# intronic read (unspliced) expression matrix
nmat <- dat$iomat;
# spanning read (intron+exon) expression matrix
smat <- dat$smat;
# 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
str(intersect(intersect(rownames(emat),rownames(nmat)),rownames(smat)))

 chr [1:5920] "Dmd" "Rbfox1" "Dlg2" "Auts2" "Ctnna2" "Rbms3" "Gpc6" "Pard3b" "Dpyd" "Hdac9" "Celf2" "Erc2" ...

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 2% of cells (by spliced expression magnitude)

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

calculating cell knn ... done
calculating convolved matrices ... done
fitting gamma coefficients ... done. succesfful fit for 8646 genes
filtered out 1335 out of 8646 genes due to low nmat-emat correlation
filtered out 968 out of 7311 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

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,deltaT=1,kCells = 5,fit.quantile = 0.02,old.fit=rvel.qf,show.gene='Chga',cell.emb=emb,cell.colors=cell.colors)

calculating convolved matrices ... done
[1] 1

Alternatively, we calculate gene-relative velocity, using k=5 cell kNN pooling, but now using entire range of expression to determine slope gamma, and using spanning reads (smat) to fit the gene offsets.

rvel <- gene.relative.velocity.estimates(emat,nmat,smat=smat,deltaT=1,kCells = 5, min.nmat.emat.slope = 0.1, min.nmat.smat.correlation = 0.1)

calculating cell knn ... done
calculating convolved matrices ... done
fitting smat-based offsets ... done
fitting gamma coefficients ... done. succesfful fit for 5920 genes
filtered out 574 out of 5920 genes due to low nmat-smat correlation
filtered out 589 out of 5346 genes due to low nmat-emat correlation
filtered out 546 out of 4757 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

Here we calculate the most basic version of velocity estimates, using relative gamma fit, without cell kNN smoothing:

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

fitting gamma coefficients ... done. succesfful fit for 8646 genes
filtered out 820 out of 8646 genes due to low nmat-emat correlation
filtered out 496 out of 7826 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

Velocity estimate based on gene structure

Genome-wide model fit:

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

filtered out 5 out of 9298 genes due to low emat levels
filtered out 1103 out of 9293 genes due to insufficient exonic or intronic lengths
filtered out 165 out of 8190 genes due to excessive nascent counts
using relative slopes for 3752 genes to fit structure-based model ... 67.3% 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 7912 genes
filtered out 930 out of 7893 genes due to low nmat-smat correlation
filtered out 911 out of 6963 genes due to low nmat-emat correlation
filtered out 388 out of 6052 genes due to low nmat-emat slope
calculating RNA velocity shift ... done
calculating extrapolated cell state ... done
re-estimated offsets for 6963 out of 8025 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

Or in 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)

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.25 seconds (sparsity = 0.904921)!
Learning embedding...
Iteration 50: error is 45.122118 (50 iterations in 1.55 seconds)
Iteration 100: error is 45.122118 (50 iterations in 1.79 seconds)
Iteration 150: error is 45.122114 (50 iterations in 1.90 seconds)
Iteration 200: error is 45.121998 (50 iterations in 1.83 seconds)
Iteration 250: error is 45.118757 (50 iterations in 1.84 seconds)
Iteration 300: error is 0.147369 (50 iterations in 1.85 seconds)
Iteration 350: error is 0.145243 (50 iterations in 1.89 seconds)
Iteration 400: error is 0.144270 (50 iterations in 1.92 seconds)
Iteration 450: error is 0.143748 (50 iterations in 1.90 seconds)
Iteration 500: error is 0.143296 (50 iterations in 1.89 seconds)
Iteration 550: error is 0.143195 (50 iterations in 1.88 seconds)
Iteration 600: error is 0.143089 (50 iterations in 2.04 seconds)
Iteration 650: error is 0.143074 (50 iterations in 1.97 seconds)
Iteration 700: error is 0.143206 (50 iterations in 1.92 seconds)
Iteration 750: error is 0.143037 (50 iterations in 1.53 seconds)
Iteration 800: error is 0.143136 (50 iterations in 1.50 seconds)
Iteration 850: error is 0.143409 (50 iterations in 1.42 seconds)
Iteration 900: error is 0.143326 (50 iterations in 1.82 seconds)
Iteration 950: error is 0.143302 (50 iterations in 2.03 seconds)
Iteration 1000: error is 0.143286 (50 iterations in 1.95 seconds)
Fitting performed in 36.42 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)

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.905338)!
Learning embedding...
Iteration 50: error is 45.060497 (50 iterations in 1.77 seconds)
Iteration 100: error is 45.060497 (50 iterations in 1.72 seconds)
Iteration 150: error is 45.060497 (50 iterations in 1.95 seconds)
Iteration 200: error is 45.060497 (50 iterations in 2.09 seconds)
Iteration 250: error is 45.060497 (50 iterations in 2.06 seconds)
Iteration 300: error is 0.176206 (50 iterations in 1.77 seconds)
Iteration 350: error is 0.161241 (50 iterations in 1.63 seconds)
Iteration 400: error is 0.162081 (50 iterations in 1.46 seconds)
Iteration 450: error is 0.161970 (50 iterations in 1.71 seconds)
Iteration 500: error is 0.161862 (50 iterations in 1.77 seconds)
Iteration 550: error is 0.162071 (50 iterations in 1.82 seconds)
Iteration 600: error is 0.162070 (50 iterations in 1.84 seconds)
Iteration 650: error is 0.161656 (50 iterations in 1.83 seconds)
Iteration 700: error is 0.161123 (50 iterations in 1.81 seconds)
Iteration 750: error is 0.162124 (50 iterations in 1.75 seconds)
Iteration 800: error is 0.161912 (50 iterations in 1.55 seconds)
Iteration 850: error is 0.162460 (50 iterations in 1.72 seconds)
Iteration 900: error is 0.161989 (50 iterations in 1.38 seconds)
Iteration 950: error is 0.161813 (50 iterations in 1.39 seconds)
Iteration 1000: error is 0.162155 (50 iterations in 1.39 seconds)
Fitting performed in 34.42 seconds.
delta norm ... done

#dev.off()

Visualization on an existing embedding

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

vel <- rvel; arrow.scale=6; 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

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 ... grid estimates ... grid.sd= 1.731182  min.arrow.size= 0.03462363  max.grid.arrow.length= 0.09156871  done

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,vel,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=500,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.