劉小澤寫于19.5.7
主要看流程，這一篇不涉及真實(shí)數(shù)據(jù)展示
總的來說，Cell Ranger主要的流程有：拆分?jǐn)?shù)據(jù) mkfastq、細(xì)胞定量 count、定量組合 aggr、調(diào)參reanalyze，還有一些小工具比如mkref、mkgtf、upload、sitecheck、mat2csv、vdj、mkvdjref、testrun
首先是mkfastq 拆分?jǐn)?shù)據(jù)
雖然這里用不到(因?yàn)槲覀兿螺d的就是fastq數(shù)據(jù))，但是為了流程的完整還是要學(xué)習(xí)一下
目的：將每個flowcell 的Illumina sequencer's base call files (BCLs)轉(zhuǎn)為fastq文件
特色： 它借鑒了Illumina出品的bcl2fastq，另外增加了：
- 將10X 樣本index名稱與四種寡核苷酸對應(yīng)起來，比如A1孔是樣本
SI-GA-A1，然后對應(yīng)的寡核苷酸是GGTTTACT, CTAAACGG, TCGGCGTC, and AACCGTAA，那么程序就會去index文件中將存在這四種寡核苷酸的fastq組合到A1這個樣本 - 提供質(zhì)控結(jié)果，包括barcode 質(zhì)量、總體測序質(zhì)量如Q30、R1和R2的Q30堿基占比、測序reads數(shù)等
- 可以使用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指的是下機(jī)的原始BCL文件目錄
# 重要的就是測序lane、樣本名稱、index等信息
samplesheet.csv文件就是illumina常規(guī)使用的，類似下面這種。它除了需要指定各種ID、name之外，還要根據(jù)不同的試劑盒版本調(diào)整[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名稱
使用簡化版的這個文件，可以識別使用的試劑盒版本，然后自行調(diào)整reads的長度信息
最后的結(jié)果就是三個文件：I1序列文件以及兩個測序文件R1、R2
目錄結(jié)構(gòu)如下：
- 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
自己分析的數(shù)據(jù)也要改成這種結(jié)構(gòu)存放，方便后續(xù)分析
小Tip--指定fastq文件位置
后續(xù)分析需要指定fastq位置，但是這些fastq文件可以由
cellranger mkfastq得到，也可以利用s Illumina'sbcl2fastq、公共數(shù)據(jù)、10X的bamtofastq,每種情況可能得到的fastq存放位置是不同的，那么如何根據(jù)不同情況進(jìn)行指定呢？
第一種情況：
利用mkfastq或者bcl2fastq生成的文件，大概長這樣
# 會有這幾種選擇方式[注意幾種參數(shù)的設(shè)置]
# 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. 所有l(wèi)anes上的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各自進(jìn)行操作
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=test_sample1,test_sample2
其實(shí)從上面的各種設(shè)置也能看出來，一開始的樣本命名規(guī)則是非常重要的
第二種情況：
也是利用mkfastq或者bcl2fastq生成的文件，但是同一個樣本的數(shù)據(jù)放在不同的目錄
# 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樣本的所有l(wèi)anes
--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樣本的所有l(wèi)anes的數(shù)據(jù)
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=MySample
# 3. 只處理Mysample樣本的lane1數(shù)據(jù)
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=MySample \
--lanes=1
第五種情況：
fastq命名方式變了，類似于這樣：
它一般是從demux流程中拆分出來的數(shù)據(jù)，但是目前被mkfastq取代，沒有好的方法，需要知道樣本相關(guān)的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數(shù)據(jù)
--fastqs=/PATH/TO/PROJECT_FOLDER \
--indices=SI-GA-A1 \
--lanes=1
# 4.利用oligo
--fastqs=/PATH/TO/PROJECT_FOLDER \
--indices=AACCGTAA,CTAAACGG,GGTTTACT,TCGGCGTC
第六種情況：
數(shù)據(jù)命名與上面完全不同，因此需要自己重命名，方式就是
# 這個在單細(xì)胞實(shí)戰(zhàn)(二)中介紹過
[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz
# 其中Read Type
# I1: Sample index read (optional)
# R1: Read 1
# R2: Read 2
分析時就可以直接調(diào)用了
# 1.所有樣本
--fastqs=/PATH/TO/PROJECT_FOLDER
# 2. 某個樣本的所有l(wèi)anes數(shù)據(jù)
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=SAMPLENAME
# 3. 某個樣本的某個lane
--sample=SAMPLENAME \
--fastqs=/PATH/TO/PROJECT_FOLDER \
--lanes=1
然后是count 細(xì)胞定量
這個過程是最重要的，它完成細(xì)胞與基因的定量，它將比對、質(zhì)控、定量都包裝了起來，內(nèi)部流程很多，但使用很簡單
先學(xué)會使用
每個版本要求的參數(shù)是不同的，尤其是V2與V3版本存在較大差異，這里先對V2進(jìn)行了解
基本上自己需要輸入的參數(shù)是：
# 這是示例，不是真實(shí)數(shù)據(jù) #
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保持一致，作為軟件識別的標(biāo)志
# expect-cells指定復(fù)現(xiàn)的細(xì)胞數(shù)量，這個要和實(shí)驗(yàn)設(shè)計(jì)結(jié)合起來
# nosecondary 只獲得表達(dá)矩陣，不進(jìn)行后續(xù)的降維、聚類和可視化分析(因?yàn)楹笃跁孕杏肦包去做)
它的輸出文件有很多
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格式數(shù)據(jù)摘要
- 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：數(shù)據(jù)分析目錄，下面又包含聚類clustering（有g(shù)raph-based & k-means）、差異分析diffexp、主成分線性降維分析pca、非線性降維tsne
- molecule_info.h5：下面進(jìn)行aggregate使用的文件
- cloupe.cloupe：官方可視化工具Loupe Cell Browser 輸入文件
一些內(nèi)置軟件和算法
基因組比對—是否在外顯子？
利用了 STAR比對工具，這款比對工具比對速度快，靈敏度高，是ENCODE、GATK推薦使用的工具，允許基因的可變剪切。比對完之后，利用GTF文件將reads溯源回外顯子區(qū)、內(nèi)含子區(qū)、基因間區(qū)：如果一條read的50%以上與外顯子有交集，那么就認(rèn)為它在外顯區(qū)；如果不在外顯子區(qū)，與內(nèi)含子有交集，那么就認(rèn)為它在內(nèi)含子區(qū)；與外顯子、內(nèi)含子都沒有交集，那么就認(rèn)為在基因間區(qū)
MAPQ 輔助判斷—在外顯子的正確率有多少？
如果reads比對到了一個外顯子區(qū)，同時也比對到了1個或多個的非外顯子區(qū)，更相信它在外顯子區(qū)，然后看MAPQ值，值越大越可信，如果MAPQ的值為255的話，那么就可以非常確定它比對到了外顯子區(qū)
MAPQ即mapping quality，告訴我們這個read比對到參考基因組上某個位置的可信度，它的公式是：
-10logP(error)，如果這個值大于30就認(rèn)為比對發(fā)生錯誤的概率是千分之一
轉(zhuǎn)錄組比對—是否特異比對？
如果上面得到的外顯子區(qū)域reads同時比對上有注釋轉(zhuǎn)錄本上的外顯子，并且在同一條鏈上，那么認(rèn)為這個reads也比對到了轉(zhuǎn)錄組；如果只比對到單個基因的注釋信息，那么認(rèn)為它是特異比對到轉(zhuǎn)錄組的(uniquely /confidently mapped )，這樣的reads才會拿來做接下來的UMI 計(jì)數(shù)
重點(diǎn)和難點(diǎn)在于自主構(gòu)建參考信息
Cell Ranger為比對和定量提供了參考基因組及注釋 pre-built human (hg19, GRCh38), mouse (mm10), and ercc92 reference packages
但是很多時候，我們需要根據(jù)自己的需要，自定義一套參考信息，但需要注意以下問題：
- 參考序列只能有很少的 overlapping gene annotations，因?yàn)閞eads比對到多個基因會導(dǎo)致流程檢測的分子數(shù)更少(它只要uniquely mapped的結(jié)果)
- 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結(jié)果中就不包含gene_biotype:pseudogene這部分
然后利用mkref構(gòu)建參考索引
# 基本使用(單個物種)
cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf
# 可以使用--nthreads使用多線程加速
# 得到的輸出結(jié)果(保存在--genome名稱的目錄中)
ls hg19
fasta/ genes/ pickle/ reference.json star/
# 如果對于多個物種組合(本文的數(shù)據(jù)其實(shí)就應(yīng)該這樣組合起來)
cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf \
--genome=mm10 --fasta=mm10.fa --genes=mm10-filtered-ensembl.gtf
# 得到的結(jié)果
ls hg19_and_mm10
fasta/ genes/ pickle/ reference.json star/
如果要增加基因信息
第一步，在fasta/genome.fa的FASTA基礎(chǔ)上增加序列信息；
第二步，在genes/genes.gtf的GTF基礎(chǔ)上增加注釋信息，注意格式
# 每一行有9列tab分隔信息
# 第一列：Chromosome 指定基因組上染色體或contig位置
# 第二列：Source 這個用處不大
# 第三列：Feature CellRanger軟件只取exon的部分
# 第四列：Start 起始位點(diǎn)(1-based)
# 第五列：End 終止位點(diǎn)(1-based)
# 第六列：Score 這個用處不大，建議用"."表示
# 第七列：Strand feature信息在基因組的+或-鏈
# 第八列：Frame 用處不大，建議“.”
# 第九列：分號分隔的鍵值對，重點(diǎn)是transcript_id 和gene_id。gene_name可選
例如：
mylocus annotation exon 100 200 . + . gene_id "mygene"; transcript_id "mygene";
第三步，使用cellranger mkref運(yùn)行更新一下
P.S. 最后得到的參考信息(包括參考基因組、注釋信息)文件結(jié)構(gòu)如下：
# 這是官網(wǎng)下載的hg38數(shù)據(jù)
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
當(dāng)處理多個生物學(xué)樣本或者一個樣本存在多個重復(fù)/文庫時，最好的操作就是先分別對每個文庫進(jìn)行單獨(dú)的count定量，然后將定量結(jié)果利用aggr組合起來
第一步 得到count結(jié)果
例如現(xiàn)在分別進(jìn)行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 ...
第二步 構(gòu)建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的路徑
第三步 運(yùn)行aggr
cellranger aggr --id=AGG123 \
--csv=AGG123_libraries.csv \
--normalize=mapped
# 結(jié)果輸出到AGG123這個目錄中
至于最后的 reanalyze ，這個屬于定制化分析，這里暫時不做探討，日后待標(biāo)準(zhǔn)化流程構(gòu)建起來，再補(bǔ)充這一部分
歡迎關(guān)注我們的公眾號～_～
我們是兩個農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com
