Analysis of single-cell RNA-seq data: building and annotating an atlas¶
This R notebook pre-processes the pbmc_1k v3 dataset from 10X Genomics with kallisto and bustools using kb, and then performs an analysis of the cell types and their marker genes.
If you use the methods in this notebook for your analysis please cite the following publications which describe the tools used in the notebook:
Melsted, P., Booeshaghi, A.S. et al. Modular and efficient pre-processing of single-cell RNA-seq. bioRxiv (2019). doi:10.1101/673285
Stuart, Butler et al. Comprehensive Integration of Single-cell Data. Cell (2019). doi:10.1016/j.cell.2019.05.031
A Python notebook implementing the same analysis is available here. See the kallistobus.tools tutorials site for additional notebooks demonstrating other analyses.
A large fraction of the running time of this notebook is in installing the Seurat R package, since it has lots of dependencies and many of them use Rcpp which results in the need to compile lots of C++ code.
'Building wheels for collected packages: loompy, numpy-groupies'
' Building wheel for loompy (setup.py): started'
' Building wheel for loompy (setup.py): finished with status \'done\''
' Created wheel for loompy: filename=loompy-3.0.6-cp36-none-any.whl size=47895 sha256=cb353db9cb57063d9db53d8ed434b4f4fdb15a77787b3eee9e60797c60c5cb57'
' Stored in directory: /root/.cache/pip/wheels/f9/a4/90/5a98ad83419732b0fba533b81a2a52ba3dbe230a936ca4cdc9'
' Building wheel for numpy-groupies (setup.py): started'
' Building wheel for numpy-groupies (setup.py): finished with status \'done\''
' Created wheel for numpy-groupies: filename=numpy_groupies-0+unknown-cp36-none-any.whl size=28042 sha256=db9efb43f40b5edece6fc869edb49644efb79cb201a2ae423a38223d51d8124a'
' Stored in directory: /root/.cache/pip/wheels/30/ac/83/64d5f9293aeaec63f9539142fc629a41af064cae1b3d8d94aa'
# Download the data from the 10x websitesystem("wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_fastqs.tar",intern=TRUE)system("tar -xvf pbmc_1k_v3_fastqs.tar",intern=TRUE)
# Slightly modified from BUSpaRse, just to avoid installing a few dependencies not used hereread_count_output<-function(dir,name){dir<-normalizePath(dir,mustWork=TRUE)m<-readMM(paste0(dir,"/",name,".mtx"))m<-Matrix::t(m)m<-as(m,"dgCMatrix")# The matrix read has cells in rowsge<-".genes.txt"genes<-readLines(file(paste0(dir,"/",name,ge)))barcodes<-readLines(file(paste0(dir,"/",name,".barcodes.txt")))colnames(m)<-barcodesrownames(m)<-genesreturn(m)}
In this plot cells are ordered by the number of UMI counts associated to them (shown on the x-axis), and the fraction of droplets with at least that number of cells is shown on the y-axis:
1
summary(tot_counts)
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 1.00 1.00 43.64 6.00 60120.00
#' Knee plot for filtering empty droplets#' #' Visualizes the inflection point to filter empty droplets. This function plots #' different datasets with a different color. Facets can be added after calling#' this function with `facet_*` functions. Will be added to the next release#' version of BUSpaRse.#' #' @param bc_rank A `DataFrame` output from `DropletUtil::barcodeRanks`.#' @return A ggplot2 object.knee_plot<-function(bc_rank){knee_plt<-tibble(rank=bc_rank[["rank"]],total=bc_rank[["total"]])%>%distinct()%>%dplyr::filter(total>0)annot<-tibble(inflection=metadata(bc_rank)[["inflection"]],rank_cutoff=max(bc_rank$rank[bc_rank$total>metadata(bc_rank)[["inflection"]]]))p<-ggplot(knee_plt,aes(total,rank))+geom_line()+geom_hline(aes(yintercept=rank_cutoff),data=annot,linetype=2)+geom_vline(aes(xintercept=inflection),data=annot,linetype=2)+scale_x_log10()+scale_y_log10()+annotation_logticks()+labs(y="Rank",x="Total UMIs")return(p)}
plot_pct_genes<-function(mat,tr2g,top_n=20,symbol="ensembl"){pct_tx<-rowSums(mat)gs<-rownames(mat)[order(-pct_tx)]df<-as.data.frame(t(mat[gs[1:20],]))df<-df%>%mutate_all(function(x)x/colSums(mat))%>%pivot_longer(everything(),names_to="gene")if (symbol=="ensembl"){df<-left_join(df,tr2g,by="gene")}else{df<-rename(df,gene_name=gene)}df%>%mutate(gene=fct_reorder(gene_name,value,.fun=median))%>%ggplot(aes(gene,value))+geom_boxplot()+labs(x="",y="Proportion of total counts")+coord_flip()}
For many barcodes, the top genes by proportion of all counts are ribosomal or mitochondrial genes. Also, the proportions plotted above seem to have some discrete values; this effect is a result of computing fractions with small denominator, which happens when droplets produce very few UMI counts.
The steps below constitute a standard analysis worklow for single-cell RNA-seq data.
1
2
# The [[ operator can add columns to object metadata. This is a great place to stash QC statspbmc[["percent.mt"]]<-PercentageFeatureSet(pbmc,pattern="^MT-")
The number of unique genes and total molecules are automatically calculated when running the CreateSeuratObject command.
The associated data is stored in the object metadata.
1
2
# Show QC metrics for the first 5 cellshead(pbmc@meta.data,5)
A data.frame: 5 × 4
orig.ident
nCount_RNA
nFeature_RNA
percent.mt
<fct>
<dbl>
<int>
<dbl>
AAACCCAAGGAGAGTA
pbmc1k
9289
3198
11.271396
AAACGCTTCAGCCCAG
pbmc1k
6483
2513
8.252352
AAAGAACAGACGACTG
pbmc1k
5011
2082
6.166434
AAAGAACCAATGGCAG
pbmc1k
3264
1555
6.893382
AAAGAACGTCTGCAAT
pbmc1k
7488
2508
6.610577
Next, we visualize some QC metrics and use the results to set filtering criteria.
1
2
3
# Visualize QC metrics as a violin plotoptions(repr.plot.width=12,repr.plot.height=6)VlnPlot(pbmc,features=c("nFeature_RNA","nCount_RNA","percent.mt"),ncol=3)
1
2
3
4
5
# FeatureScatter is typically used to visualize feature-feature relationships, but can be used# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.plot1<-FeatureScatter(pbmc,feature1="nCount_RNA",feature2="percent.mt")plot2<-FeatureScatter(pbmc,feature1="nCount_RNA",feature2="nFeature_RNA")CombinePlots(plots=list(plot1,plot2))
After removing unwanted cells from the dataset, the next step is to normalize the data. A standard choice is LogNormalize which normalizes the UMI counts for each cell by the total counts, multiplies this by a scale factor (10,000 by default), and finally log-transforms the result. Normalized values are stored in pbmc[["RNA"]]@data.
We recommend the preprint
- Breda, J., Zavolan, M. and van Nimwegen, E. Bayesian inference of the gene expression states of single cells from scRNA-seq data. bioRxiv (2019). doi.org/10.1101/2019.12.28.889956
For clarity, in this previous line of code (and in future commands), we provide the default values for certain parameters in the function call. However, this isn’t required and the same behavior can be achieved with:
To identify a subset of genes that exhibit high cell-to-cell variation in the dataset we apply a procedure implemented in the FindVariableFeatures function. By default, it returns 2,000 genes per dataset. These will be used in downstream analysis.
Seurat documentation describes the method used to find highly variable genes here as such:
First, fits a line to the relationship of log(variance) and log(mean) using local polynomial regression (loess). Then standardizes the feature values using the observed mean and expected variance (given by the fitted line). Feature variance is then calculated on the standardized values after clipping to a maximum (see clip.max parameter).
options(repr.plot.width=9,repr.plot.height=6)pbmc<-FindVariableFeatures(pbmc,selection.method="vst",nfeatures=2000)# Identify the 10 most highly variable genestop10<-head(VariableFeatures(pbmc),10)# plot variable features with and without labelsplot1<-VariableFeaturePlot(pbmc,log=FALSE)LabelPoints(plot=plot1,points=top10,repel=TRUE)
Warning message:“Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.Please use `as_label()` or `as_name()` instead.[90mThis warning is displayed once per session.[39m”When using repel, set xnudge and ynudge to 0 for optimal results
Next, we apply a linear transformation (‘scaling’) that is a standard pre-processing step prior to dimensional reduction techniques like PCA. The ScaleData function shifts the expression of each gene, so that the mean expression across cells is 0 and the variance across cells is 1
This step gives equal weight to genes in downstream analyses, so that highly-expressed genes do not dominate. The results of this are stored in pbmc[["RNA"]]@scale.data
We apply this only to the genes identified as highly variable:
1
# pbmc <- ScaleData(pbmc)
The scaling does not affect PCA or clustering results. However, Seurat heatmaps (produced as shown below with DoHeatmap) require genes in the heatmap to be scaled so that highly-expressed genes don’t dominate. To make sure we don’t leave any genes out of the heatmap later, we are scaling all genes in this tutorial.
In Seurat v2 we also use the ScaleData function to remove unwanted sources of variation from a single-cell dataset. For example, we could ‘regress out’ heterogeneity associated with (for example) cell cycle stage, or mitochondrial contamination. These features are still supported in ScaleData in Seurat v3, i.e.:
To overcome the extensive technical noise in any single feature for scRNA-seq data, one can cluster cells based on their PCA projections, with each PC essentially representing a ‘metafeature’ that combines information across a correlated feature set.
A common heuristic method generates an ‘Elbow plot’: a ranking of principle components based on the percentage of variance explained by each one (ElbowPlot function). In this example, we can observe an ‘elbow’ around PC9-10, suggesting that the majority of true signal is captured in the first 10 PCs.
We cluster cells using the Louvain algorithm (a default in Seurat), which iteratively group cells together, with the goal of optimizing the standard modularity function. The FindClusters function implements this procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.4-1.2 typically returns good results for single-cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters can be found using the Idents function.
tSNE and UMAP can be used to visualize and explore non-linear aspects of high-dimensional data. Here we apply these methods to the PC projection of the data (with same dimension as used for clustering).
UMAP (UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction) is a manifold learning technique that can also be used to visualize cells. It was published in:
McInnes, Leland, John Healy, and James Melville. "Umap: Uniform manifold approximation and projection for dimension reduction." arXiv preprint arXiv:1802.03426 (2018).
t-SNE is a non-linear dimensionality reduction technique described in:
Maaten, Laurens van der, and Geoffrey Hinton. "Visualizing data using t-SNE." Journal of machine learning research 9.Nov (2008): 2579-2605.
1
pbmc<-RunUMAP(pbmc,dims=1:10,verbose=FALSE)
Warning message:“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metricTo use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'This message will be shown once per session”
1
2
3
# note that you can set `label = TRUE` or use the LabelClusters function to help label# individual clustersDimPlot(pbmc,reduction="umap")
Finding differentially expressed features (cluster biomarkers)¶
A key follow-up step to clustering cells is to find gene markers that are associated with them. We used Seurat's FindAllMarkers function which automates the process for all clusters.
The min.pct argument requires a feature to be detected at a minimum percentage in either of the two groups of cells, and the thresh.test argument requires a feature to be differentially expressed (on average) by some amount between the two groups. You can set both of these to 0, but with a dramatic increase in time - since this will test a large number of features that are unlikely to be highly discriminatory. As another option to speed up these computations, max.cells.per.ident can be set. This will downsample each identity class to have no more cells than whatever this is set to. While there is generally going to be a loss in power, the speed increases can be significiant and the most highly differentially expressed features will likely still rise to the top.
Several methods for differential expression are supported by Seurat. The default is Wilcoxon rank sum test.
1
2
3
# find markers for every cluster compared to all remaining cells, report only the positive onespbmc.markers<-FindAllMarkers(pbmc,test.use="wilcox",only.pos=TRUE,min.pct=0.25,logfc.threshold=0.25)
Seurat includes several tools for visualizing marker expression. VlnPlot (shows expression probability distributions across clusters), and FeaturePlot (visualizes feature expression on a tSNE or PCA plot) are our most commonly used visualizations. We also suggest exploring RidgePlot, CellScatter, and DotPlot as additional methods to view your dataset.