公共數(shù)據(jù)庫(kù)中的SRA 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)究竟包含了哪些數(shù)據(jù)讲婚？CellRanger怎樣利用10x平臺(tái)下機(jī)數(shù)據(jù)進(jìn)行下游一系列分析笙蒙？這篇文章簡(jiǎn)單記錄CellRanger 包括的主要分析步驟凯旋，純理論。
SRA 原始數(shù)據(jù)轉(zhuǎn)fastq
公共數(shù)據(jù)庫(kù)的SRA 數(shù)據(jù)需要借助fastq-dump 轉(zhuǎn)為fastq文件样漆，然后進(jìn)行質(zhì)控晌缘、CellRanger定量等操作。相較于普通轉(zhuǎn)錄組數(shù)據(jù)亿鲜，原始SRA數(shù)據(jù)會(huì)得到3個(gè)fastq文件允蜈，分別是Library 的Index(8bp)文件，包括長(zhǎng)度為26bp 的Barcode（16bp）和UMI(10bp)的Read1文件蒿柳，和測(cè)序reads文件饶套。
conda install -c bioconda sra-tools ## 安裝軟件
wkd=/home/project/single-cell/MCC
cd $wkd/raw/P2586-4
cat SRR_Acc_List-2586-4.txt |while read i
do
time fastq-dump --gzip --split-files -A $i ${i}.sra && echo "** ${i}.sra to fastq done **"
done
### 單細(xì)胞數(shù)據(jù)參數(shù)為 --split-files 而不是 --split-3
i7 sample index (library barcode）
是加到Illumina測(cè)序接頭上的，保證多個(gè)測(cè)序文庫(kù)可以在同一個(gè)flow-cell上或者同一個(gè)lane上進(jìn)行混合測(cè)序(multiplexed)垒探。它的作用就是在CellRanger的mkfastq 功能中體現(xiàn)出來的妓蛮，它自動(dòng)識(shí)別樣本index名稱（例如：SA-GA-A1），將具有相同4種oligo的fq文件組合在一起圾叼，表示同一個(gè)樣本蛤克。它保證了一個(gè)測(cè)序lane上可以容納多個(gè)樣本
Barcode
是10X特有的，用來區(qū)分GEMs夷蚊，也就是對(duì)細(xì)胞做了一個(gè)標(biāo)記构挤。一般在拆分混樣測(cè)序數(shù)據(jù)（demultiplexing）這個(gè)過程后進(jìn)行操作，當(dāng)然這也很符合原文的操作惕鼓。
UMI
UMI就是Unique Molecular Identifier筋现，由4-10個(gè)隨機(jī)核苷酸組成，在mRNA反轉(zhuǎn)錄后箱歧，進(jìn)入到文庫(kù)中矾飞，每一個(gè)mRNA隨機(jī)連上一個(gè)UMI，根據(jù)PCR結(jié)果可以計(jì)數(shù)不同的UMI呀邢，最終統(tǒng)計(jì)mRNA的數(shù)量凰慈。它的主要作用是，處理PCR 擴(kuò)增偏差驼鹅，因?yàn)槠鹗嘉膸?kù)很小時(shí)需要的PCR擴(kuò)增次數(shù)就越多微谓，因?yàn)樵饺菀滓霐U(kuò)增誤差。
fastq文件更名
為什么要更名输钩？CellRanger 定量過程輸入文件指定命名格式豺型。
如何更名？下圖格式：
# 比如买乃，將原來的SRR7692286_1.fastq.gz改成SRR7692286_S1_L001_I1_001.fastq.gz
# 依次類推姻氨，將原來_2的改成R1，將_3改成R2
cat SRR_Acc_List-9245-3.txt | while read i ;do (mv ${i}_1*.gz ${i}_S1_L001_I1_001.fastq.gz;mv ${i}_2*.gz ${i}_S1_L001_R1_001.fastq.gz;mv ${i}_3*.gz ${i}_S1_L001_R2_001.fastq.gz);done
fastq 文件質(zhì)控
主要查看Q30比例剪验。
# 以P2586-4為例
mkdir -p $wkd/qc
cd $wkd/qc
find $wkd/raw/P2586-4 -name '*R1*.gz'>P2586-4-id-1.txt
find $wkd/raw/P2586-4 -name '*R2*.gz'>P2586-4-id-2.txt
cat P2586-4-id-1.txt P2586-4-id-2.txt >P2586-4-id-all.txt
cat P2586-4-id-all.txt| xargs fastqc -t 20 -o ./
CellRanger 流程
它主要包括四個(gè)主要基因表達(dá)分析流程：
-
cellranger mkfastq ： 它借鑒了Illumina的
bcl2fastq肴焊，可以將一個(gè)或多個(gè)lane中的混樣測(cè)序樣本按照index標(biāo)簽生成樣本對(duì)應(yīng)的fastq文件前联。即對(duì)下機(jī)數(shù)據(jù)base calling files轉(zhuǎn)為fastq文件。 -
cellranger count ：利用
mkfastq生成的fq文件娶眷，進(jìn)行比對(duì)(基于STAR)似嗤、過濾、UMI計(jì)數(shù)届宠。利用細(xì)胞的barcode生成gene-barcode矩陣烁落，然后進(jìn)行樣本分群、基因表達(dá)分析 -
cellranger aggr ：接受
cellranger count的輸出數(shù)據(jù)豌注，將同一組的不同測(cè)序樣本的表達(dá)矩陣整合在一起伤塌，比如tumor組原來有4個(gè)樣本，PBMC組有兩個(gè)樣本轧铁，現(xiàn)在可以使用aggr生成最后的tumor和PBMC兩個(gè)矩陣每聪，并且進(jìn)行標(biāo)準(zhǔn)化去掉測(cè)序深度的影響 -
cellranger reanalyze ：接受
cellranger count或cellranger aggr生成的gene-barcode矩陣，使用不同的參數(shù)進(jìn)行降維齿风、聚類药薯。屬于定制化分析。
它的結(jié)果主要是包含有細(xì)胞信息的BAM, MEX, CSV, HDF5 and HTML文件
CellRanger 軟件安裝及測(cè)試實(shí)戰(zhàn)操作參考：CellRanger走起(四) CellRanger流程概覽聂宾。包括安裝詳細(xì)指南及其他一些小技巧果善。
拆分?jǐn)?shù)據(jù) mkfastq诊笤、細(xì)胞定量 count系谐、定量組合 aggr
mkfastq 拆分
目的：將每個(gè)flowcell 的Illumina sequencer's base call files (BCLs)轉(zhuǎn)為fastq文件
特色： 它借鑒了Illumina出品的bcl2fastq，另外增加了：
- 將10X 樣本index名稱與四種寡核苷酸對(duì)應(yīng)起來讨跟，比如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等信息
參考文章CellRanger走起(四) CellRanger流程概覽 中解釋了幾種不同方式輸出fastq文件的不同存放目錄向叉，可根據(jù)實(shí)際分析自行查找。
count 定量
這個(gè)過程是最重要的嗦董，它完成細(xì)胞與基因的定量母谎，它將比對(duì)、質(zhì)控京革、定量都包裝了起來奇唤，內(nèi)部流程很多幸斥，但使用很簡(jiǎn)單。每個(gè)版本要求的參數(shù)是不同的咬扇，尤其是V2與V3版本存在較大差異甲葬，這里先對(duì)V2進(jìn)行了解。
# 這是示例冗栗，不是真實(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é)合起來
# 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!
輸出文件從上到下依次來看：
- web_summary.html：官方說明 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：過濾掉的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 輸入文件
count 定量的重點(diǎn)和難點(diǎn)在于參考基因組索引的構(gòu)建烦粒，可詳細(xì)參考CellRanger走起(四) CellRanger流程概覽休溶，文章中還介紹了怎樣自行構(gòu)建參考基因組索引，及一比對(duì)過程的重要細(xì)節(jié)信息扰她。
aggr
當(dāng)處理多個(gè)生物學(xué)樣本或者一個(gè)樣本存在多個(gè)重復(fù)/文庫(kù)時(shí)兽掰，最好的操作就是先分別對(duì)每個(gè)文庫(kù)進(jìn)行單獨(dú)的count定量，然后將定量結(jié)果利用aggr組合起來.
- 得到count結(jié)果
- 構(gòu)建Aggregation CSV
- 運(yùn)行aggr
### step1
$ 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 ...
## step2
# 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的路徑
### step3
cellranger aggr --id=AGG123 \
--csv=AGG123_libraries.csv \
--normalize=mapped
# 結(jié)果輸出到AGG123這個(gè)目錄中
cellranger count的結(jié)果
count結(jié)果一般放在out目錄下徒役，主要有summary和analysis兩大類
Summary
是一個(gè)HTML格式文件孽尽，包括實(shí)驗(yàn)捕獲的細(xì)胞數(shù)目、檢測(cè)到的基因數(shù)目忧勿、測(cè)序數(shù)據(jù)的產(chǎn)出與質(zhì)量統(tǒng)計(jì)杉女、參考基因組的比對(duì)情況。
幾個(gè)指標(biāo)可以關(guān)注一下：
左上部分中鸳吸，包括了reads數(shù)熏挎、barcodes數(shù)、UMI层释、index婆瓜、Q30等統(tǒng)計(jì)值
左下是reads比對(duì)的比例，包括基因間區(qū)、外顯子廉白、內(nèi)含子个初，如果比對(duì)率太低(一般認(rèn)為外顯子的比對(duì)率要在60%以上)
右上圖是利用barcodes上的UMI標(biāo)簽分布來估計(jì)細(xì)胞數(shù)，綠/藍(lán)色表示細(xì)胞猴蹂，灰色表示背景院溺，其中Y軸表示每個(gè)barcode對(duì)應(yīng)的UMI數(shù)量，X軸是一定數(shù)量的UMI序列所對(duì)應(yīng)barcode的數(shù)量磅轻，比如上圖中有1000個(gè)barcodes含有10k個(gè)UMI珍逸，細(xì)胞過濾就是通過這個(gè)圖來展示的。
首先明確聋溜，barcode用來區(qū)分細(xì)胞谆膳，UMI用來區(qū)分轉(zhuǎn)錄本；其次撮躁，barcodes數(shù)量時(shí)要大于細(xì)胞數(shù)量的（以保證每個(gè)細(xì)胞都會(huì)有barcode來進(jìn)行區(qū)分）如果圖線陡降說明
另參考文章CellRanger走起(二) CellRanger out 結(jié)果 中還記錄了一個(gè)手動(dòng)構(gòu)建物種GTF文件的案例嘗試漱病，及相應(yīng)的檢驗(yàn)方法，可以根據(jù)需要自行學(xué)習(xí)把曼。
總之杨帽，無論是從公共數(shù)據(jù)庫(kù)下載單細(xì)胞SRA數(shù)據(jù)，還是來自10x 平臺(tái)的原始下機(jī)數(shù)據(jù)嗤军，或者先經(jīng)過fastq-dump轉(zhuǎn)為fastq文件注盈，然后質(zhì)控，然后CellRanger count 定量叙赚，或者首先自行拆分?jǐn)?shù)據(jù)mkfastq老客，然后構(gòu)建索引進(jìn)行細(xì)胞定量，及可選的定量組合等纠俭。CellRanger的總體流程見本文記錄沿量，但是其中的各個(gè)細(xì)節(jié)仍需要自行實(shí)踐總結(jié)浪慌。
作者：新欣enjoy
鏈接：http://www.reibang.com/p/5157ab9f6977
