前面说了cellranger安装和一些简单的使用,我们继续看一下cellranger的结果。
一文学会单细胞转录组的CellRanger(一) - 简书 (jianshu.com)
一、HTML报告解读
cellranger count 输出一个名为web_summary.html
的交互式 HTML 文件,其中包含一些汇总的指标和二次分析结果。如果在运行期间检测到问题,则此页面上会显示warning或error。Cell Ranger 故障排除文档中提供了有关警报的详细信息,可在文档中进行查询。
注:如设置了--nosecondary参数则没有二次分析的结果
1、基本功能
图是HTML的部分截图(1)可以通过单击 HTML 文件中的“Sequencing”、“Mapping”和“Cells”旁边的?图标,来查看指标的含义,点击后出现如下图的信息:
(2)可以通过单击 HTML 文件左上角选项卡中的Summary来查看一些汇总指标。汇总指标描述了测序质量和检测到的细胞的各种特征。点击Analysis查看二次分析结果。
2、重要指标含义
(1)、下图显示的是检测到的估计细胞数、每个细胞的平均reads数和每个细胞检测到的基因中位数。(2)、“Cells"部分的barcode-rank-plot的分布图。y 轴是映射到每个barcode的 UMI 计数值,x 轴是低于该值的barcode数。急剧下降表示与cell关联的barcode与背景区关联的barcode之间具有良好的分离。由于barcode可以根据其UMI计数或RNA谱与细胞相关联,因此图形的某些区域可以同时包含细胞相关和背景相关barcode。图形的颜色表示与cell关联的barcode的局部密度。
(3)、测序饱和度Sequencing Saturation: 每个样本过滤后细胞的reads数占总reads数(含背景)的百分比,反应测序数据利用率,当值达到80%以上,代表测序reads覆盖了所有mRNA。
(4)、若Reads Mapped Confidently to Intronic Regions的值高于40%,是因为包含pre-mRNA。
(5)、Fraction Reads in Cells :理想情况高于70%,数据质量则比较好。若值较低,可能是以下两种情况导致。第一,细胞裂解或死亡,cell ranger将其判断为背景RNA,导致样品中背景RNA含量过高;第二,大量具有低 RNA 含量的细胞,算法并未将其判断为cell导致的。 后一种情况可以通过检查数据来确定适当的细胞计数并使用 --force-cells 来解决。
二、filtered_feature_bc_matrix和raw_feature_bc_matrix文件夹
两个文件夹下的内容是一样的,是稀疏矩阵,包含三个文件。文件的区别是数量上的,顾名思义,一个是原始的数据,一个是过滤后的结果。后面的分析一般使用filtered_feature_bc_matrix文件夹下的矩阵。
filtered_feature_bc_matrix
├── barcodes.tsv.gz
├── features.tsv.gz
└── matrix.mtx.gz
1、 barcodes.tsv.gz文件存放的细胞的barcodes标签,可以理解为每个barcode标签代表一个cell。
$ less barcodes.tsv.gz|head
AAACCCAAGACGCATG-1
AAACCCACACAAGCCC-1
AAACCCACACTTACAG-1
AAACCCAGTAACTTCG-1
AAACCCAGTCTGTTAG-1
AAACCCAGTGAGATAT-1
AAACCCAGTTCGTACA-1
AAACCCATCCGAGATT-1
AAACGAAAGACTTGTC-1
AAACGAAAGCTCACTA-1
2、 features.tsv.gz文件存放的gene id和对应的gene name。
$ less features.tsv.gz|head
ENSG00000243485 MIR1302-2HG Gene Expression
ENSG00000237613 FAM138A Gene Expression
ENSG00000186092 OR4F5 Gene Expression
ENSG00000238009 AL627309.1 Gene Expression
ENSG00000239945 AL627309.3 Gene Expression
ENSG00000239906 AL627309.2 Gene Expression
ENSG00000241860 AL627309.5 Gene Expression
ENSG00000241599 AL627309.4 Gene Expression
ENSG00000286448 AP006222.2 Gene Expression
ENSG00000236601 AL732372.1 Gene Expression
3、matrix.mtx.gz文件存放的对应barcode和features的对应count值。前两行是一些版本信息等。第三行依次是features数量总和,barcode数量总和,count值总和,第四行开始是对应的结果。比如 31 1 4 对应的是在features文件第31行的基因,1对应的barcode文件中的第1行信息,4对应检测到的count值。
$ less matrix.mtx.gz|head
%%MatrixMarket matrix coordinate integer general
%metadata_json: {"software_version": "cellranger-4.0.0", "format_version": 2}
36601 9620 11561775
31 1 4
74 1 2
81 1 1
84 1 1
87 1 1
114 1 1
171 1 3
三、analysis文件夹
软件默认会进行二次分析,结果中会出现analysis文件。如果后续用其他软件代替这部分内容(比如Seurat、Scanpy等),可以设置了--nosecondary参数,省略二次分析。以下是analysis文件夹的结构:
clustering
├── graphclust
├── kmeans_10_clusters
├── kmeans_2_clusters
├── kmeans_3_clusters
├── kmeans_4_clusters
├── kmeans_5_clusters
├── kmeans_6_clusters
├── kmeans_7_clusters
├── kmeans_8_clusters
└── kmeans_9_clusters
diffexp
├── graphclust
├── kmeans_10_clusters
├── kmeans_2_clusters
├── kmeans_3_clusters
├── kmeans_4_clusters
├── kmeans_5_clusters
├── kmeans_6_clusters
├── kmeans_7_clusters
├── kmeans_8_clusters
└── kmeans_9_clusters
pca
└── 10_components
tsne
└── 2_components
umap
└── 2_components
1、pca 降维
在对细胞进行聚类之前,在归一化过滤的矩阵上运行主成分分析 (PCA),以减少features(基因)维度的数量。仅将基因表达特征用作 PCA 特征。PCA 分析生成五个输出文件。第一个是每个cell对前 N 个主分量的投影。默认情况下 N=10。
$ head -2 analysis/pca/gene_expression_10_components/projection.csv
Barcode,PC-1,PC-2,PC-3,PC-4,PC-5,PC-6,PC-7,PC-8,PC-9,PC-10
AAACAAGCACCATACT-1,18.55496347631502,-8.428877305709332,3.7717969735420835,-0.61215157678172,-1.0987614379684771,2.194733668965279,-2.6595895212967386,-2.8703699622639114,1.867229094193604,0.2658532968798859
第二个文件是一个分量矩阵,它指示每个特征对每个主分量的贡献(荷载)。未包含在 PCA 分析中的要素的所有载荷值都设置为零。
$ head -2 analysis/pca/gene_expression_10_components/components.csv
PC,ENSG00000228327,ENSG00000237491,ENSG00000177757,ENSG00000225880,...,ENSG00000160310
1,-0.0044,0.0039,-0.0024,-0.0016,...,-0.0104
第三个文件包含选择用于主成分计算的具有最高离散的要素的gene id。
$ head -5 analysis/pca/gene_expression_10_components/features_selected.csv
Feature
1,ENSG00000167723
2,ENSG00000179029
3,ENSG00000196544
4,ENSG00000141499
第四个文件记录每个主成分解释的总方差比例。 在选择重要的主成分数量时,查看很有用,当数字变化平缓时, 后续 PC 在数据中的意义不大。
$ head -5 analysis/pca/gene_expression_10_components/variance.csv
PC,Proportion.Variance.Explained
1,0.0056404970744118104
2,0.0038897311237809061
3,0.0028803714818085419
4,0.0020830581822081206
最后一个文件列出了每个要素,按平均表达式对要素进行分箱后的归一化离散程度,用于度量每个特征的可变性。
$ head -5 analysis/pca/gene_expression_10_components/dispersion.csv
Feature,Normalized.Dispersion
ENSG00000228327,2.0138970131886671
ENSG00000237491,1.3773662040549017
ENSG00000177757,-0.28102027567224191
ENSG00000225880,1.9887312950109921
2、t-SNE
运行 PCA 后,运行 t-distributed Stochastic Neighbor Embedding(t-SNE) 将数据在一个2D的维度进行可视化。
$ head -5 analysis/tsne/gene_expression_2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300
3、UMAP
运行 PCA 后,运行Uniform Manifold Approximation and Projection(UMAP)将数据在一个2D的维度进行可视化。
$ head -5 analysis/umap/gene_expression_2_components/projection.csv
Barcode,UMAP-1,UMAP-2
AAACCTGAGAATAGGG-1,0.5974335,1.320372
AAACCTGAGAGCTGGT-1,2.2277818,-0.52756095
AAACCTGAGCGTTGCC-1,2.675832,1.1010709
AAACCTGCACGGACAA-1,2.7049212,-3.1494563
4、clustering 聚类
运行聚类分析,根据具有相似表达谱的细胞在 PCA 空间中的投影,将它们分组在一起。cellranger使用了两中方法:
- Graph-based
图聚类算法包括两步:首先用PCA降维的数据构建一个细胞间的k近邻稀疏矩阵,即将一个细胞与其欧式距离上最近的k个细胞聚为一类,然后在此基础上用Louvain算法进行模块优化,旨在找到图中高度连接的模块。最后通过层次聚类将位于同一区域内没有差异表达基因(B-H adjusted p-value 低于0.05)的cluster进一步融合,重复该过程直到没有clusters可以合并。因为它不需要预先指定数量的聚类,只需要运行一次。 - K-Means
k-means算法随机在PCA降维的空间中适当选取k个聚类质心点,对于每一个细胞计算其应该属于的cluster,然后对于每一个cluster重新计算该cluster的质心,重复该过程直到收敛。注意这里K-means针对 K=2,...,N 的许多值运行,其中 K 对应于聚类数。默认情况下 N=10(与图聚类算法的k意义不同),质心代表对属于同一个cluster的细胞中心点的猜测。k-means可说是最简单、最经典的聚类算法。
$ ls analysis/clustering
gene_expression_graphclust
gene_expression_kmeans_10_clusters
gene_expression_kmeans_2_clusters
gene_expression_kmeans_3_clusters
gene_expression_kmeans_4_clusters
gene_expression_kmeans_5_clusters
gene_expression_kmeans_6_clusters
gene_expression_kmeans_7_clusters
gene_expression_kmeans_8_clusters
gene_expression_kmeans_9_clusters
5、diffexp差异表达
cellranger还会生成一个表,指示每个聚类中相对于所有其他聚类中哪些要素以差异方式表示。对于每个特征和每个聚类 i,我们计算三个值:
- The mean expression of this feature in cluster i (i.e., across cells assigned to cluster i)
- The log2 fold-change of this feature's mean expression in cluster i relative to all other cells
- A p-value denoting significance of this feature's expression in cluster i relative to cells in other clusters. P-values within each cluster are adjusted for false discovery rate to account for the number of hypotheses (i.e., number of features) being tested.
$ head -5 analysis/diffexp/gene_expression_kmeans_3_clusters/differential_expression.csv
Feature ID,Feature Name,Cluster 1 Mean UMI Counts,Cluster 1 Log2 fold change,Cluster 1 Adjusted p value,Cluster 2 Mean UMI Counts,Cluster 2 Log2 fold change,Cluster 2 Adjusted p value,Cluster 3 Mean UMI Counts,Cluster 3 Log2 fold change,Cluster 3 Adjusted p value
ENSG00000228327,RP11-206L10.2,0.0056858989363338264,2.6207666981569986,0.00052155805898912184,0.0,-0.75299726644507814,0.64066099091888962,0.00071455453829430329,-2.3725403666493312,0.0043023680184636837
ENSG00000237491,RP11-206L10.9,0.00012635330969630726,-0.31783275717885928,0.40959138980118809,0.0,3.8319652342760779,0.11986963938734894,0.0,0.56605908868652577,0.39910771338768203
ENSG00000177757,FAM87B,0.0,-2.9027952579000154,0.0,0.0,3.2470027335549219,0.19129034227967889,0.00071455453829430329,3.1510215894076818,0.0
ENSG00000225880,LINC00115,0.0003790599290889218,-5.71015017995762,8.4751637615375386e-28,0.20790015775229512,7.965820981010868,1.3374521290889345e-46,0.0017863863457357582,-2.2065304152104019,0.00059189960914085744