单细胞实战(四) Cell Ranger流程概览 / 开普饭

实战演练

理论知识学再好，能付诸实践灵活运用才行，所以我们常强调知行合一，实践出真知。实战演练这个栏目就是带大家从头到尾完整复现单细胞文献分析流程。好了，干货多，屁话少，我们来看实战流程。

希望大家能有所收获！

前情回顾

总的来说，Cell Ranger主要的流程有：拆分数据 mkfastq、细胞定量 count、定量组合 aggr、调参reanalyze，还有一些小工具比如mkref、mkgtf、upload、sitecheck、mat2csv、vdj、mkvdjref、testrun

首先是mkfastq拆分数据

虽然这里用不到(因为我们下载的就是fastq数据)，但是为了流程的完整还是要学习一下

目的：将每个flowcell 的Illumina sequencer's base call files (BCLs)转为fastq文件

特色： 它借鉴了Illumina出品的bcl2fastq，另外增加了：

将10X 样本index名称与四种寡核苷酸对应起来，比如A1孔是样本SI-GA-A1，然后对应的寡核苷酸是GGTTTACT, CTAAACGG, TCGGCGTC, and AACCGTAA ，那么程序就会去index文件中将存在这四种寡核苷酸的fastq组合到A1这个样本
提供质控结果，包括barcode 质量、总体测序质量如Q30、R1和R2的Q30碱基占比、测序reads数等
可以使用10X简化版的样本信息表

它的示意流程：

两种使用方式：

# 第一种
$ cellranger mkfastq --id=bcl \
--run=/path/to/bcl \
--samplesheet=samplesheet-1.2.0.csv
# 第二种
$ cellranger mkfastq --id=bcl \
--run=/path/to/bcl \
--csv=simple-1.2.0.csv
# 其中id指定输出目录的名称，run指的是下机的原始BCL文件目录
# 重要的就是测序lane、样本名称、index等信息

samplesheet.csv文件就是illumina常规使用的，类似下面这种。它除了需要指定各种ID、name之外，还要根据不同的试剂盒版本调整[Reads]长度

V2试剂盒R1序列长度为26bp(包括16bp的barcode+10bp的UMI)，R2为98bp；V3试剂盒R1序列长度为28bp(包括16bp的barcode+12bp的UMI)，R2为91bp

还有一种10X定制的简单化的csv文件，例如：

Lane,Sample,Index
1,test_sample,SI-GA-A3
# 其中第一列指定lane ID，第二列是样本名称，第三列是index名称

使用简化版的这个文件，可以识别使用的试剂盒版本，然后自行调整reads的长度信息

最后的结果就是三个文件：I1序列文件以及两个测序文件R1、R2

目录结构如下：

- tiny-bcl/outs/fastq_path/bcl/
- Sample1
- Sample1_S1_L001_I1_001.fastq.gz
- Sample1_S1_L001_R1_001.fastq.gz
- Sample1_S1_L001_R2_001.fastq.gz

自己分析的数据也要改成这种结构存放，方便后续分析

小Tip--指定fastq文件位置

后续分析需要指定fastq位置，但是这些fastq文件可以由cellranger mkfastq得到，也可以利用s Illumina's bcl2fastq 、公共数据、10X的bamtofastq ,每种情况可能得到的fastq存放位置是不同的，那么如何根据不同情况进行指定呢？

第一种情况：

利用mkfastq或者bcl2fastq生成的文件，大概长这样

# 会有这几种选择方式[注意几种参数的设置]
# 1.所有mkfastq生成的样本
--fastqs=MKFASTQ_ID/outs/fastq_path
# 2. 多个flowcell生成的所有样本
--fastqs=MKFASTQ_ID/outs/fastq_path1,MKFASTQ_ID/outs/fastq_path2
# 3.所有bcl2fastq 生成的样本
--fastqs=/PATH/TO/bcl2fastq_output
# 4. 所有lanes上的test_sample1样本
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample1
# 5. lane1上的test_sample1样本
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample1 \
--lanes=1
# 6. 将test_sample1和test_sample2各自进行操作
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample1,test_sample2

其实从上面的各种设置也能看出来，一开始的样本命名规则是非常重要的

第二种情况：

也是利用mkfastq或者bcl2fastq生成的文件，但是同一个样本的数据放在不同的目录

# 1. 将所有SI-GA-A1样本的reads组合
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=SI-GA-A1_1,SI-GA-A1_2,SI-GA-A1_3,SI-GA-A1_4
# 2. 只处理SI-GA-A1样本的第一个index样本
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=SI-GA-A1_1

第三种情况：

也是利用mkfastq或者bcl2fastq生成的文件，但和Reports、Stats在同一个目录

# 1. mkfastq得到的所有样本
--fastqs=MKFASTQ_ID/outs/fastq_path
# 2. bcl2fastq得到的所有样本
--fastqs=/PATH/TO/bcl2fastq_output
# 3. test_sample样本的所有lanes
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample
# 4. test_sample样本的lane1
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample \
--lanes=1

第四种情况：

使用 mkfastq or bcl2fastq 得到的fastq文件和Report、Stats不在同一个目录，但命名方式与之前一样，这个目录中只能看到fastq文件

# 1.处理所有样本
--fastqs=/PATH/TO/PROJECT_FOLDER
# 2. 处理Mysample样本的所有lanes的数据
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=MySample
# 3. 只处理Mysample样本的lane1数据
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=MySample \
--lanes=1

第五种情况：

fastq命名方式变了，类似于这样：

它一般是从demux流程中拆分出来的数据，但是目前被mkfastq取代，没有好的方法，需要知道样本相关的index或者oligos

# 1.所有样本
--fastqs=/PATH/TO/PROJECT_FOLDER
# 2.所有SI-GA-A1样本
--fastqs=/PATH/TO/PROJECT_FOLDER \
--indices=SI-GA-A1
# 3.所有SI-GA-A1样本的lane1数据
--fastqs=/PATH/TO/PROJECT_FOLDER \
--indices=SI-GA-A1 \
--lanes=1
# 4.利用oligo
--fastqs=/PATH/TO/PROJECT_FOLDER \
--indices=AACCGTAA,CTAAACGG,GGTTTACT,TCGGCGTC

第六种情况：

数据命名与上面完全不同，因此需要自己重命名，方式就是

# 这个在单细胞实战(二)中介绍过
[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz
# 其中Read Type
# I1: Sample index read (optional)
# R1: Read 1
# R2: Read 2

分析时就可以直接调用了

# 1.所有样本
--fastqs=/PATH/TO/PROJECT_FOLDER
# 2. 某个样本的所有lanes数据
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=SAMPLENAME
# 3. 某个样本的某个lane
--sample=SAMPLENAME \
--fastqs=/PATH/TO/PROJECT_FOLDER \
--lanes=1

然后是count细胞定量

这个过程是最重要的，它完成细胞与基因的定量，它将比对、质控、定量都包装了起来，内部流程很多，但使用很简单

先学会使用

每个版本要求的参数是不同的，尤其是V2与V3版本存在较大差异，这里先对V2进行了解

基本上自己需要输入的参数是：

# 这是示例，不是真实数据 #
cellranger count --id=sample345 \
--transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0 \
--fastqs=/home/scRNA/runs/HAWT7ADXX/outs/fastq_path \
--sample=mysample \
--expect-cells=1000 \
--nosecondary
# id指定输出文件存放目录名
# transcriptome指定与CellRanger兼容的参考基因组
# fastqs指定mkfastq或者自定义的测序文件
# sample要和fastq文件的前缀中的sample保持一致，作为软件识别的标志
# expect-cells指定复现的细胞数量，这个要和实验设计结合起来
# nosecondary 只获得表达矩阵，不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做)

它的输出文件有很多

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

Pipestance completed successfully!

从上到下依次来看：

web_summary.html：官方说明 summary HTML file
metrics_summary.csv：CSV格式数据摘要
possorted_genome_bam.bam：比对文件
possorted_genome_bam.bam.bai：索引文件
filtered_gene_bc_matrices：是重要的一个目录，下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件
filtered_feature_bc_matrix.h5：过滤掉的barcode信息HDF5 format
raw_feature_bc_matrix：原始barcode信息
raw_feature_bc_matrix.h5：原始barcode信息HDF5 format
analysis：数据分析目录，下面又包含聚类clustering（有graph-based & k-means）、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne
molecule_info.h5：下面进行aggregate使用的文件
cloupe.cloupe：官方可视化工具Loupe Cell Browser 输入文件

一些内置软件和算法

基因组比对—是否在外显子？

利用了 STAR比对工具，这款比对工具比对速度快，灵敏度高，是ENCODE、GATK推荐使用的工具，允许基因的可变剪切。比对完之后，利用GTF文件将reads溯源回外显子区、内含子区、基因间区：如果一条read的50%以上与外显子有交集，那么就认为它在外显区；如果不在外显子区，与内含子有交集，那么就认为它在内含子区；与外显子、内含子都没有交集，那么就认为在基因间区

MAPQ 辅助判断—在外显子的正确率有多少？

如果reads比对到了一个外显子区，同时也比对到了1个或多个的非外显子区，更相信它在外显子区，然后看MAPQ值，值越大越可信，如果MAPQ的值为255的话，那么就可以非常确定它比对到了外显子区

MAPQ即mapping quality，告诉我们这个read比对到参考基因组上某个位置的可信度，它的公式是：-10logP(error)，如果这个值大于30就认为比对发生错误的概率是千分之一

转录组比对—是否特异比对？

如果上面得到的外显子区域reads同时比对上有注释转录本上的外显子，并且在同一条链上，那么认为这个reads也比对到了转录组；如果只比对到单个基因的注释信息，那么认为它是特异比对到转录组的(uniquely /confidently mapped )，这样的reads才会拿来做接下来的UMI 计数

重点和难点在于自主构建参考信息

Cell Ranger为比对和定量提供了参考基因组及注释 pre-built human (hg19, GRCh38), mouse (mm10), and ercc92 reference packages

但是很多时候，我们需要根据自己的需要，自定义一套参考信息，但需要注意以下问题：

参考序列只能有很少的 overlapping gene annotations，因为reads比对到多个基因会导致流程检测的分子数更少(它只要uniquely mapped的结果)
FASTA与GTF比对和STAR兼容，GTF文件的第三列（feature type）必须有exon

首先利用mkgtf过滤GTF文件

先从 ENSEMBL或UCSC上下载，然后使用mkgtf

cellranger mkgtf input.gtf output.gtf --attribute=key:allowable_value
# 其中键值对可以指定多个，比如
$ cellranger mkgtf Homo_sapiens.GRCh38.ensembl.gtf Homo_sapiens.GRCh38.ensembl.filtered.gtf \
--attribute=gene_biotype:protein_coding \
--attribute=gene_biotype:lincRNA \
--attribute=gene_biotype:antisense \
--attribute=gene_biotype:IG_LV_gene \
--attribute=gene_biotype:IG_V_gene \
--attribute=gene_biotype:IG_V_pseudogene \
--attribute=gene_biotype:IG_D_gene \
--attribute=gene_biotype:IG_J_gene \
--attribute=gene_biotype:IG_J_pseudogene \
--attribute=gene_biotype:IG_C_gene \
--attribute=gene_biotype:IG_C_pseudogene \
--attribute=gene_biotype:TR_V_gene \
--attribute=gene_biotype:TR_V_pseudogene \
--attribute=gene_biotype:TR_D_gene \
--attribute=gene_biotype:TR_J_gene \
--attribute=gene_biotype:TR_J_pseudogene \
--attribute=gene_biotype:TR_C_gene
# 这样得到的Homo_sapiens.GRCh38.ensembl.filtered.gtf结果中就不包含gene_biotype:pseudogene这部分

然后利用mkref构建参考索引

# 基本使用(单个物种)
cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf
# 可以使用--nthreads使用多线程加速
# 得到的输出结果(保存在--genome名称的目录中)
ls hg19
fasta/ genes/ pickle/ reference.json star/

# 如果对于多个物种组合(本文的数据其实就应该这样组合起来)
cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf \
--genome=mm10 --fasta=mm10.fa --genes=mm10-filtered-ensembl.gtf
# 得到的结果
ls hg19_and_mm10
fasta/ genes/ pickle/ reference.json star/

如果要增加基因信息

参考链接：https://kb.10xgenomics.com/hc/en-us/articles/115003327112-How-can-we-add-genes-to-a-reference-package-for-Cell-Ranger-

第一步，在fasta/genome.fa的FASTA基础上增加序列信息；

第二步，在genes/genes.gtf的GTF基础上增加注释信息，注意格式

# 每一行有9列tab分隔信息
# 第一列：Chromosome 指定基因组上染色体或contig位置
# 第二列：Source 这个用处不大
# 第三列：Feature CellRanger软件只取exon的部分
# 第四列：Start 起始位点(1-based)
# 第五列：End 终止位点(1-based)
# 第六列：Score 这个用处不大，建议用"."表示
# 第七列：Strand feature信息在基因组的+或-链
# 第八列：Frame 用处不大，建议“.”
# 第九列：分号分隔的键值对，重点是transcript_id 和gene_id。gene_name可选
例如：
mylocus annotation exon 100 200 . + . gene_id "mygene"; transcript_id "mygene";

第三步，使用cellranger mkref运行更新一下

P.S. 最后得到的参考信息(包括参考基因组、注释信息)文件结构如下：

4 directories, 21 files

多个文库的整合 aggr

当处理多个生物学样本或者一个样本存在多个重复/文库时，最好的操作就是先分别对每个文库进行单独的count定量，然后将定量结果利用aggr组合起来

第一步得到count结果

例如现在分别进行3个定量流程

$ cellranger count --id=LV123 ...
... wait for pipeline to finish ...
$ cellranger count --id=LB456 ...
... wait for pipeline to finish ...
$ cellranger count --id=LP789 ...
... wait for pipeline to finish ...

第二步构建Aggregation CSV

就像这样：

# AGG123_libraries.csv
library_id,molecule_h5
LV123,/opt/runs/LV123/outs/molecule_info.h5
LB456,/opt/runs/LB456/outs/molecule_info.h5
LP789,/opt/runs/LP789/outs/molecule_info.h5
# 其中
# molecule_h5：文件molecule_info.h5 file的路径

第三步运行aggr

cellranger aggr --id=AGG123 \
--csv=AGG123_libraries.csv \
--normalize=mapped
# 结果输出到AGG123这个目录中

至于最后的 reanalyze ，这个属于定制化分析，这里暂时不做探讨，日后待标准化流程构建起来，再补充这一部分

单细胞实战(四) Cell Ranger流程概览