cellranger使用的初步探索（2）理解cellranger count输出文件

在上一篇笔记的结尾处，提到了cellranger count这一步，之前几天我一直是用的服务器里128G和200G的内存试着运行这一步，并且限制了cores的数量（64），但是在24小时之内都没有处理完SRR7722937这一个样品的fastq文件（在服务器里跑程序需要设置运行时间）。这次我取消了内存以及cores的数量限制（这时cellranger会耗尽几乎所有的可用内存，大概将近300个G），并且把运行时间的限制从24小时改成了5天：

$ cellranger count \
          --id=Tumor937 
          --transcriptome=/gpfs/home/practice/10_genomics_genome/GRCh38 \ 
          --fastqs=/gpfs/home/practice \ 
          --sample=SRR7722937
          --expect-cells=10000

这次的尝试成功的跑完了一个样品，总共运行时间为2天8个小时。打开log文件可以看到程序运行过程中记录了很多条记录，大概有上百行记录，在记录的最后几行，会出现运行的总结，说明你的cellranger count运行成功：

Outputs:
- Run summary HTML:                      /gpfs/home/practice/Tumor937/outs/web_summary.html
- Run summary CSV:                       /gpfs/home/practice/Tumor937/outs/metrics_summary.csv
- BAM:                                   /gpfs/home/practice/Tumor937/outs/possorted_genome_bam.bam
- BAM index:                             /gpfs/home/practice/Tumor937/outs/possorted_genome_bam.bam.bai
- Filtered gene-barcode matrices MEX:    /gpfs/home/practice/Tumor937/outs/filtered_gene_bc_matrices
- Filtered gene-barcode matrices HDF5:   /gpfs/home/practice/Tumor937/outs/filtered_gene_bc_matrices_h5.h5
- Unfiltered gene-barcode matrices MEX:  /gpfs/home/practice/Tumor937/outs/raw_gene_bc_matrices
- Unfiltered gene-barcode matrices HDF5: /gpfs/home/practice/Tumor937/outs/raw_gene_bc_matrices_h5.h5
- Secondary analysis output CSV:         /gpfs/home/practice/Tumor937/outs/analysis
- Per-molecule read information:         /gpfs/home/practice/Tumor937/outs/molecule_info.h5
- Loupe Cell Browser file:               /gpfs/home/practice/Tumor937/outs/cloupe.cloupe

Waiting 6 seconds for UI to do final refresh.
Pipestance completed successfully!

Saving pipestance info to Tumor937/Tumor937.mri.tgz

打开我们在cellranger count代码里设置的文件夹Tumor937，里面有以下文件：

$ ll
total 4820
-rw------- 1 fy04 fy04     267 Aug  4 20:35 _cmdline
-rw------- 1 fy04 fy04   48206 Aug  4 20:35 _filelist
-rw------- 1 fy04 fy04  725168 Aug  4 20:35 _finalstate
-rw------- 1 fy04 fy04     655 Aug  2 12:01 _invocation
-rw------- 1 fy04 fy04       5 Aug  2 12:01 _jobmode
-rw------- 1 fy04 fy04  246834 Aug  4 20:35 _log
-rw------- 1 fy04 fy04   46606 Aug  2 12:01 _mrosource
drwx------ 5 fy04 fy04    8192 Aug  4 20:35 outs #重点关注这个文件夹
-rw------- 1 fy04 fy04  324837 Aug  4 20:35 _perf
drwx------ 5 fy04 fy04    8192 Aug  4 20:34 SC_RNA_COUNTER_CS
-rw------- 1 fy04 fy04   12160 Aug  4 20:35 _sitecheck
-rw------- 1 fy04 fy04       2 Aug  2 12:01 _tags
-rw------- 1 fy04 fy04      51 Aug  4 20:35 _timestamp
-rw------- 1 fy04 fy04 3164754 Aug  4 20:35 Tumor937.mri.tgz
-rw------- 1 fy04 fy04      36 Aug  2 12:01 _uuid
-rw------- 1 fy04 fy04  102558 Aug  4 20:35 _vdrkill
-rw------- 1 fy04 fy04      61 Aug  2 12:01 _versions

再打开“outs”文件夹：

$ ll
total 3264857
drwx------ 6 fy04 fy04       8192 Aug  4 20:35 analysis
-rw------- 1 fy04 fy04   27898707 Aug  4 20:35 cloupe.cloupe #这个文件可以使用Loupe Cell Browser软件打开
drwx------ 3 fy04 fy04       8192 Aug  4 20:35 filtered_gene_bc_matrices
-rw------- 1 fy04 fy04    7247168 Aug  4 16:04 filtered_gene_bc_matrices_h5.h5 #过滤后的基因-barcode矩阵，HDF5格式
-rw------- 1 fy04 fy04        684 Aug  4 20:34 metrics_summary.csv #csv格式的结果总结
-rw------- 1 fy04 fy04  105930133 Aug  4 17:47 molecule_info.h5 #之后如果需要aggr整合几个单细胞测序的样品，这个文件是aggr步骤的输入文件
-rw------- 1 fy04 fy04 3182359671 Aug  4 15:59 possorted_genome_bam.bam #比对到基因组和转录本上的reads，并且带有barcode注释信息
-rw------- 1 fy04 fy04    3949040 Aug  4 16:00 possorted_genome_bam.bam.bai #possorted_genome_bam.bam的索引文件
drwx------ 3 fy04 fy04       8192 Aug  4 20:35 raw_gene_bc_matrices #没有过滤的基因-barcode矩阵，包括所有的barcode
-rw------- 1 fy04 fy04   12085113 Aug  4 16:02 raw_gene_bc_matrices_h5.h5
-rw------- 1 fy04 fy04    3530757 Aug  4 20:34 web_summary.html #网页版的运行结果总结，可以使用任何一个浏览器查看

解读cellranger count的输出文件

（1）web_summary.html文件

一旦cellranger count运行完毕，你可以通过这个文件在浏览器里查看你的结果总结。你也可以使用Loupe Cell Browser软件查看“outs”文件夹里的.cloupe文件。下面这个截屏就是总结，里面记录了测序质量和检测到的细胞的一些特征值。包括检测到的细胞数量，每个细胞的平均reads数，每个细胞检测到的平均基因数（上方的绿色的三个数字）。右面的曲线图显示的是barcode的计数分布，以及哪些barcode对应的相关细胞。Y轴是map到每一个barcode的UMI的计数数值，X轴是对应到这个数值的barcode的数量。这条曲线要有一个急剧下降的趋势才是好的结果，说明细胞barcode在细胞和“空样品”间的差异非常明显。这个曲线图下面的“Fraction reads in cells”这一项比例需要大于70%，说明数据才是质量好的。

上面这个网页的左上角，有一个“ANALYSIS”的选项，点开看，会出现下面这个页面：

这个页面里包括以下几个部分：
（1）t-SNE降维分析
（2）自动生成的聚类分析（根据表达谱里相似表达的基因）
（3）一个差异基因列表，是根据自动生成的聚类来生成的差异基因
（4）一个测序饱和度的曲线
（5）每个细胞里测得的基因数量的曲线

上面图里左边的图是一个二维的t-SNE，每一个点就是一个细胞，根据细胞里含有的UMI数量来区分颜色。这张图也说明了细胞内RNA含量，也通常与细胞大小有关。红色的点代表这个细胞有着更多的RNA含量。在这个坐标系里，两个离得很近的点有着更为相似的基因表达模式（相比于两个离得较远的点）。

而上面的右边的是一个聚类的图，你可以通过改变右上角的clustering type，来改变聚类的参数：

这个页面的中部的表就是差异基因列表：

再往下两张图就是测序饱和度（通常也指示文库的复杂度），以及平均每个细胞里测得的基因，这两张图都可以说明测序深度：

（2）bam文件

在“outs”文件夹里，有一个bam文件，是根据我们的fastq文件生成的：

$ samtools view possorted_genome_bam.bam | head -5
SRR7722937.377377       16      1       10232   255     4S94M   *       0       0       TTTTCCCTATCCCTAACCCTAACCCTATCCCCTTACCCCTAACCCTAACCCTAACCCTAACCCCAACCCCAACCCCAACCCCAACCCCTACCCCAACC    .<.<..<<...<.<...G<...A.<....AGG<.GGGGA<.<GGGG<<GGGGGGGIGGG<GGA<.<AAA.<<GA<A<<AGAAA<GA<<..<.<..<AG    NH:i:1  HI:i:1  AS:i:72 nM:i:10 RE:A:I  BC:Z:TTTCATGA   QT:Z:GGGGGIII CR:Z:CTTACCGAGTGTACCT   CY:Z:GGGGGIIIIIIIIIII   CB:Z:CTTACCGAGTGTACCT-1 UR:Z:TTTCCATCGT UY:Z:IIIIIIIIIG UB:Z:TTTCCATCGT       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.10911097     256     1       11201   0       1S97M   *       0       0       GGCTTGCTCACGGTGCTGTGCCAGGGCGCCCCCTGCTGGCGACTAGGGCAACTGCAGGGCTCTCTTGCTTAGACTGGTGGCCAGCGCCCCCTGCTGGC    .AG...GAGGIIGGGGGGGGIGGGGIIIIGGGGGIIIIIIIIGGGIGGIGGGIIIIIGIAGAAGAAGGGGGGGGIGGGGGGG<<.AAAGGIGAGAAA<    NH:i:5  HI:i:2  AS:i:93 nM:i:1  RE:A:I  BC:Z:ACGTCCCT   QT:Z:GGGGGIII CR:Z:CACCACTGTCCATGAT   CY:Z:GGGGGIIIIGIIIIII   CB:Z:CACCACTGTCCATGAT-1 UR:Z:GAGCGCGGGC UY:Z:IIIIIIIIII UB:Z:GAGCGCGGGC       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.27657393     256     1       11208   0       98M     *       0       0       CACGGTGCTGTGCCAGGGCGCCCCCTGCTGGCGACTAGGGCAACTGCAGGGCTCTCTTGCTTAGAGTGGTGGCCAGCGCCCCCTGCTGGCGCCGGGGC    <.A...7A..<<.<GGG.G<AAAGGGIIIGGIIIGIGAA<GG.GGAGGGGGGGGIGGIIIGGAGG.<GGGGAGAAGAGGI<AAAGGAGGGAAA.<A..    NH:i:5  HI:i:3  AS:i:96 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:AGTGAGGGTCAGGACA   CY:Z:GGGGGIIIIIIIIIII   CB:Z:AGTGAGGGTCAGGACA-1 UR:Z:CAGAACACCA UY:Z:IIIIIIIIII UB:Z:CAGAACACCA       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.33049902     256     1       11285   0       98M     *       0       0       GCCCCCTGCTGGCGCCGGGGCACTGCAGGGCCCTCTTGCTTACTGTATAGTGGTGGCACGCCGCCTGCTGGCAGCTAGGGACATTGCAGGGTCCTCTT    AGAGAA<AG.<AA.AGGGGGIGGGGGGGGGGGGGIGGIIGIIGGG<<AGGGAGGGGGAGGGAG<AGGGGGGIGG.GAGAAAA.GGGG.AGGGGGAA.G    NH:i:5  HI:i:2  AS:i:96 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:ATTACTCTCCCTTGTG   CY:Z:GGGGAIIIIIIIIIII   CB:Z:ATTACTCTCCCTTGTG-1 UR:Z:CAATATCCCC UY:Z:IGGGIIIIII UB:Z:CAATATCCCC       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.33235779     256     1       11292   0       97M1S   *       0       0       GCTGGCGCCGGGGCACTGCAGGGCCCTCTTGCTTACTGTATAGTGGTGGCACGCCGCCTGCTGGCAGCTAGGGACATTGCAGGGTCCTCTTGCTCAAA    GGGA...<A<<AGGGGGGAGAGGAGGIIIGIIIIGGGIGGIIIIIIIIGIGGIGGGGAGGGGGGIGGAGGA<GGGGGGGGGGI<GGGGGIIIAAGGG.    NH:i:6  HI:i:3  AS:i:95 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:TGCCCTAGTAGTAGTA   CY:Z:GGGGAIIIIIGGGIII   CB:Z:TGCCCTAGTAGTAGTA-1 UR:Z:GTCAGAAGTG UY:Z:IIIIIIIIII UB:Z:GTCAGAAGTG       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1

这个bam文件有25列，前11列是标准的bam文件的内容，后面的14列是cellranger count生成的bam文件特有的，这25列分别是：

其中第10列是Read序列，第11列是read的测序质量；第19列是barcode序列（CR），第20列是barcode测序质量（CY）；第22列是UMI序列（UR），第23列是UMI测序质量（UY）

（3）矩阵文件

运行cellranger count后会自动生成的二级分析结果，但是一般来说，我们通常是拿到基因表达矩阵后，自己在R里进行后续的分析。
参考文章：https://davetang.org/muse/2018/08/09/getting-started-with-cell-ranger/

在“outs”文件夹里，有一个名为raw_gene_bc_matrices的文件夹：

#包含3个文件：
$ tree
.
└── GRCh38
    ├── barcodes.tsv
    ├── genes.tsv
    └── matrix.mtx

这3个文件你可以在R里用DropletUtils包直接加载，生成一个SingleCellExperiment对象：

> BiocManager::install("DropletUtils")
> library("DropletUtils")
> sce <- read10xCounts('raw_gene_bc_matrices/GRCh38/')
> sce
class: SingleCellExperiment 
dim: 33694 737280 
metadata(1): Samples
assays(1): counts
rownames(33694): ENSG00000243485 ENSG00000237613 ...
  ENSG00000277475 ENSG00000268674
rowData names(2): ID Symbol
colnames: NULL
colData names(2): Sample Barcode
reducedDimNames(0):
spikeNames(0):
altExpNames(0):

barcodeRanks函数可以把barcode按照总UMI计数来排序：

> br.out <- barcodeRanks(counts(sce))
> library("dplyr")
> br.out.df <- as.data.frame(br.out)
> br.out.df$barcode <- colData(sce)$Barcode
> br.out.df$knee <- br.out@metadata[["knee"]]
> br.out.df$inflection <- br.out@metadata[["inflection"]]
> br.out.df %>%
  filter(rank <= 10) %>%
  arrange(rank)

    rank total   fitted            barcode knee inflection
1     1 40819       NA CAGGTGCAGCGCTTAT-1 2843       1282
2     2 39401       NA AACTCAGGTTACGGAG-1 2843       1282
3     3 38740       NA CTGATAGTCACCAGGC-1 2843       1282
4     4 36129       NA TTAACTCAGACACTAA-1 2843       1282
5     5 34162       NA GTATCTTCATGCTAGT-1 2843       1282
6     6 33527       NA TTCTCAATCGTACGGC-1 2843       1282
7     7 32814       NA CTCGTACAGCACGCCT-1 2843       1282
8     8 31075       NA ATGTGTGGTCTGGTCG-1 2843       1282
9     9 30447 30753.25 CTCGAGGTCGCGATCG-1 2843       1282
10   10 30439 29994.15 AACCGCGCAGACGCCT-1 2843       1282

我们可以自己画一个barcode vs UMI的点图，就像网页版结果总结里的曲线图那样：

> library(ggplot2)
> x_knee <- br.out.df %>% filter(br.out.df$total > br.out.df$knee) %>% arrange(total) %>% select(rank) %>% head(1)
> x_inflection <- br.out.df %>% filter(br.out.df$total > br.out.df$inflection) %>% > arrange(total) %>% select(rank) %>% head(1)
> padding <- length(br.out$rank) / 10

> p1 <- ggplot(br.out.df, aes(x = rank, y = total)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  theme_bw() +
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 14),
        title = element_text(size = 16)) +
  geom_hline(yintercept = br.out$knee, linetype = 2, colour = "dodgerblue") +
  geom_hline(yintercept = br.out$inflection, linetype = 2, colour = "forestgreen") +
  labs(x = "Rank", y = "Total", title = "Barcode Rank vs Total UMI") +
  annotate("text", label = paste0("Knee (", x_knee, ")"), x = x_knee$rank + padding, y = br.out.df$knee, size = 5) +
  annotate("text", label = paste0("Inflection (", x_inflection, ")"), x = x_inflection$rank + padding, y = br.out.df$inflection, size = 5)
> p1

NOTE：我运行的cellranger版本是2.2的，现在最新版已经是cellranger3了，所以结果部分可能会有些出入。