劉小澤寫(xiě)于19.5.7
主要看流程，這一篇不涉及真實(shí)數(shù)據(jù)展示
總的來(lái)說(shuō)，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í)一下
目的：將每個(gè)flowcell 的Illumina sequencer's base call files (BCLs)轉(zhuǎn)為fastq文件
特色： 它借鑒了Illumina出品的bcl2fastq慎冤，另外增加了：
- 將10X 樣本index名稱與四種寡核苷酸對(duì)應(yīng)起來(lái)疼燥，比如A1孔是樣本
SI-GA-A1，然后對(duì)應(yīng)的寡核苷酸是GGTTTACT, CTAAACGG, TCGGCGTC, and AACCGTAA粪薛，那么程序就會(huì)去index文件中將存在這四種寡核苷酸的fastq組合到A1這個(gè)樣本 - 提供質(zhì)控結(jié)果悴了，包括barcode 質(zhì)量、總體測(cè)序質(zhì)量如Q30违寿、R1和R2的Q30堿基占比、測(cè)序reads數(shù)等
- 可以使用10X簡(jiǎn)化版的樣本信息表
它的示意流程：
兩種使用方式：
# 第一種
$ 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文件目錄
# 重要的就是測(cè)序lane藤巢、樣本名稱、index等信息
samplesheet.csv文件就是illumina常規(guī)使用的息罗，類似下面這種掂咒。它除了需要指定各種ID、name之外，還要根據(jù)不同的試劑盒版本調(diào)整[Reads]長(zhǎng)度
V2試劑盒R1序列長(zhǎng)度為26bp(包括16bp的barcode+10bp的UMI)绍刮，R2為98bp温圆；
V3試劑盒R1序列長(zhǎng)度為28bp(包括16bp的barcode+12bp的UMI)，R2為91bp
還有一種10X定制的簡(jiǎn)單化的csv文件孩革，例如：
Lane,Sample,Index
1,test_sample,SI-GA-A3
# 其中第一列指定lane ID岁歉，第二列是樣本名稱，第三列是index名稱
使用簡(jiǎn)化版的這個(gè)文件膝蜈，可以識(shí)別使用的試劑盒版本锅移，然后自行調(diào)整reads的長(zhǎng)度信息
最后的結(jié)果就是三個(gè)文件：I1序列文件以及兩個(gè)測(cè)序文件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生成的文件肺素，大概長(zhǎng)這樣
# 會(huì)有這幾種選擇方式[注意幾種參數(shù)的設(shè)置]
# 1.所有mkfastq生成的樣本
--fastqs=MKFASTQ_ID/outs/fastq_path
# 2. 多個(gè)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è)置也能看出來(lái)，一開(kāi)始的樣本命名規(guī)則是非常重要的
第二種情況：
也是利用mkfastq或者bcl2fastq生成的文件深浮，但是同一個(gè)樣本的數(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樣本的第一個(gè)index樣本
--fastqs=MKFASTQ_ID/outs/fastq_path \
--sample=SI-GA-A1_1
第三種情況：
也是利用mkfastq或者bcl2fastq生成的文件压怠，但和Reports、Stats在同一個(gè)目錄
# 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不在同一個(gè)目錄菌瘫，但命名方式與之前一樣，這個(gè)目錄中只能看到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流程中拆分出來(lái)的數(shù)據(jù)雨让，但是目前被mkfastq取代，沒(méi)有好的方法忿等，需要知道樣本相關(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ù)命名與上面完全不同栖忠，因此需要自己重命名，方式就是
# 這個(gè)在單細(xì)胞實(shí)戰(zhàn)(二)中介紹過(guò)
[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz
# 其中Read Type
# I1: Sample index read (optional)
# R1: Read 1
# R2: Read 2
分析時(shí)就可以直接調(diào)用了
# 1.所有樣本
--fastqs=/PATH/TO/PROJECT_FOLDER
# 2. 某個(gè)樣本的所有l(wèi)anes數(shù)據(jù)
--fastqs=/PATH/TO/PROJECT_FOLDER \
--sample=SAMPLENAME
# 3. 某個(gè)樣本的某個(gè)lane
--sample=SAMPLENAME \
--fastqs=/PATH/TO/PROJECT_FOLDER \
--lanes=1
然后是count 細(xì)胞定量
這個(gè)過(guò)程是最重要的贸街，它完成細(xì)胞與基因的定量庵寞，它將比對(duì)、質(zhì)控薛匪、定量都包裝了起來(lái)捐川，內(nèi)部流程很多，但使用很簡(jiǎn)單
先學(xué)會(huì)使用
每個(gè)版本要求的參數(shù)是不同的逸尖，尤其是V2與V3版本存在較大差異古沥，這里先對(duì)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或者自定義的測(cè)序文件
# sample要和fastq文件的前綴中的sample保持一致，作為軟件識(shí)別的標(biāo)志
# expect-cells指定復(fù)現(xiàn)的細(xì)胞數(shù)量岩齿，這個(gè)要和實(shí)驗(yàn)設(shè)計(jì)結(jié)合起來(lái)
# nosecondary 只獲得表達(dá)矩陣太颤，不進(jìn)行后續(xù)的降維、聚類和可視化分析(因?yàn)楹笃跁?huì)自行用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!
從上到下依次來(lái)看：
- web_summary.html：官方說(shuō)明 summary HTML file
- metrics_summary.csv：CSV格式數(shù)據(jù)摘要
- possorted_genome_bam.bam：比對(duì)文件
- possorted_genome_bam.bam.bai：索引文件
- filtered_gene_bc_matrices：是重要的一個(gè)目錄盹沈，下面又包含了 barcodes.tsv.gz龄章、features.tsv.gz、matrix.mtx.gz襟诸，是下游Seurat瓦堵、Scater、Monocle等分析的輸入文件
- filtered_feature_bc_matrix.h5：過(guò)濾掉的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)置軟件和算法
基因組比對(duì)—是否在外顯子惋鸥？
利用了 STAR比對(duì)工具，這款比對(duì)工具比對(duì)速度快悍缠，靈敏度高卦绣，是ENCODE、GATK推薦使用的工具飞蚓，允許基因的可變剪切滤港。比對(duì)完之后，利用GTF文件將reads溯源回外顯子區(qū)趴拧、內(nèi)含子區(qū)溅漾、基因間區(qū)：如果一條read的50%以上與外顯子有交集，那么就認(rèn)為它在外顯區(qū)著榴；如果不在外顯子區(qū)添履，與內(nèi)含子有交集，那么就認(rèn)為它在內(nèi)含子區(qū)脑又；與外顯子暮胧、內(nèi)含子都沒(méi)有交集，那么就認(rèn)為在基因間區(qū)
MAPQ 輔助判斷—在外顯子的正確率有多少问麸？
如果reads比對(duì)到了一個(gè)外顯子區(qū)往衷，同時(shí)也比對(duì)到了1個(gè)或多個(gè)的非外顯子區(qū)，更相信它在外顯子區(qū)严卖，然后看MAPQ值炼绘，值越大越可信，如果MAPQ的值為255的話妄田，那么就可以非常確定它比對(duì)到了外顯子區(qū)
MAPQ即mapping quality俺亮，告訴我們這個(gè)read比對(duì)到參考基因組上某個(gè)位置的可信度，它的公式是：
-10logP(error)疟呐，如果這個(gè)值大于30就認(rèn)為比對(duì)發(fā)生錯(cuò)誤的概率是千分之一
轉(zhuǎn)錄組比對(duì)—是否特異比對(duì)脚曾？
如果上面得到的外顯子區(qū)域reads同時(shí)比對(duì)上有注釋轉(zhuǎn)錄本上的外顯子，并且在同一條鏈上启具，那么認(rèn)為這個(gè)reads也比對(duì)到了轉(zhuǎn)錄組本讥；如果只比對(duì)到單個(gè)基因的注釋信息，那么認(rèn)為它是特異比對(duì)到轉(zhuǎn)錄組的(uniquely /confidently mapped )鲁冯，這樣的reads才會(huì)拿來(lái)做接下來(lái)的UMI 計(jì)數(shù)
重點(diǎn)和難點(diǎn)在于自主構(gòu)建參考信息
Cell Ranger為比對(duì)和定量提供了參考基因組及注釋 pre-built human (hg19, GRCh38), mouse (mm10), and ercc92 reference packages
但是很多時(shí)候拷沸，我們需要根據(jù)自己的需要，自定義一套參考信息薯演，但需要注意以下問(wèn)題：
- 參考序列只能有很少的 overlapping gene annotations撞芍，因?yàn)閞eads比對(duì)到多個(gè)基因會(huì)導(dǎo)致流程檢測(cè)的分子數(shù)更少(它只要uniquely mapped的結(jié)果)
- FASTA與GTF比對(duì)和STAR兼容，GTF文件的第三列（feature type）必須有exon跨扮，過(guò)濾后的GTF只包含有注釋的基因類型
首先利用mkgtf過(guò)濾GTF文件
先從 ENSEMBL或UCSC上下載序无，然后使用mkgtf
cellranger mkgtf input.gtf output.gtf --attribute=key:allowable_value
# 其中鍵值對(duì)可以指定多個(gè)，比如
$ 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)建參考索引
# 基本使用(單個(gè)物種)
cellranger mkref --genome=hg19 --fasta=hg19.fa --genes=hg19-filtered-ensembl.gtf
# 可以使用--nthreads使用多線程加速
# 得到的輸出結(jié)果(保存在--genome名稱的目錄中)
ls hg19
fasta/ genes/ pickle/ reference.json star/
# 如果對(duì)于多個(gè)物種組合(本文的數(shù)據(jù)其實(shí)就應(yīng)該這樣組合起來(lái))
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 這個(gè)用處不大
# 第三列：Feature CellRanger軟件只取exon的部分
# 第四列：Start 起始位點(diǎn)(1-based)
# 第五列：End 終止位點(diǎn)(1-based)
# 第六列：Score 這個(gè)用處不大哟玷，建議用"."表示
# 第七列：Strand feature信息在基因組的+或-鏈
# 第八列：Frame 用處不大，建議“.”
# 第九列：分號(hào)分隔的鍵值對(duì)一也，重點(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
多個(gè)文庫(kù)的整合 aggr
當(dāng)處理多個(gè)生物學(xué)樣本或者一個(gè)樣本存在多個(gè)重復(fù)/文庫(kù)時(shí)讼渊，最好的操作就是先分別對(duì)每個(gè)文庫(kù)進(jìn)行單獨(dú)的count定量，然后將定量結(jié)果利用aggr組合起來(lái)
第一步 得到count結(jié)果
例如現(xiàn)在分別進(jìn)行3個(gè)定量流程
$ 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這個(gè)目錄中
至于最后的 reanalyze 尊剔，這個(gè)屬于定制化分析爪幻，這里暫時(shí)不做探討，日后待標(biāo)準(zhǔn)化流程構(gòu)建起來(lái)须误，再補(bǔ)充這一部分
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩挨稿，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)京痢、通俗易懂的生信知識(shí)平臺(tái)奶甘。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com
