总的来说,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'sbcl2fastq
、公共数据、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,过滤后的GTF只包含有注释的基因类型
首先利用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. 最后得到的参考信息(包括参考基因组、注释信息)文件结构如下:
# 这是官网下载的hg38数据 refdata-cellranger-GRCh38-1.2.0 |-- [ 222] README.BEFORE.MODIFYING |-- [4.0K] fasta | `-- [2.9G] genome.fa |-- [4.0K] genes | `-- [1.3G] genes.gtf |-- [4.0K] pickle | `-- [ 58M] genes.pickle |-- [ 424] reference.json |-- [4.0K] star | |-- [3.0G] Genome | |-- [8.0G] SA | |-- [1.5G] SAindex | |-- [1.2K] chrLength.txt | |-- [1.9K] chrName.txt | |-- [3.0K] chrNameLength.txt | |-- [2.1K] chrStart.txt | |-- [ 37M] exonGeTrInfo.tab | |-- [ 15M] exonInfo.tab | |-- [526K] geneInfo.tab | |-- [ 909] genomeParameters.txt | |-- [9.1M] sjdbInfo.txt | |-- [7.1M] sjdbList.fromGTF.out.tab | |-- [7.1M] sjdbList.out.tab | `-- [9.4M] transcriptInfo.tab `-- [ 6] version 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这个目录中