單細胞測序分析上游cellranger軟件入門介紹

公共數(shù)據(jù)庫中的SRA 單細胞轉(zhuǎn)錄組數(shù)據(jù)究竟包含了哪些數(shù)據(jù)？CellRanger怎樣利用10x平臺下機數(shù)據(jù)進行下游一系列分析？這篇文章簡單記錄CellRanger 包括的主要分析步驟，純理論。

SRA 原始數(shù)據(jù)轉(zhuǎn)fastq

公共數(shù)據(jù)庫的SRA 數(shù)據(jù)需要借助fastq-dump 轉(zhuǎn)為fastq文件，然后進行質(zhì)控、CellRanger定量等操作。相較于普通轉(zhuǎn)錄組數(shù)據(jù)，原始SRA數(shù)據(jù)會得到3個fastq文件，分別是Library 的Index(8bp)文件，包括長度為26bp 的Barcode（16bp）和UMI(10bp)的Read1文件，和測序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
### 單細胞數(shù)據(jù)參數(shù)為 --split-files 而不是 --split-3
image
image

i7 sample index (library barcode）
是加到Illumina測序接頭上的，保證多個測序文庫可以在同一個flow-cell上或者同一個lane上進行混合測序(multiplexed)。它的作用就是在CellRanger的mkfastq 功能中體現(xiàn)出來的，它自動識別樣本index名稱（例如：SA-GA-A1），將具有相同4種oligo的fq文件組合在一起，表示同一個樣本。它保證了一個測序lane上可以容納多個樣本

Barcode
是10X特有的，用來區(qū)分GEMs，也就是對細胞做了一個標(biāo)記。一般在拆分混樣測序數(shù)據(jù)（demultiplexing）這個過程后進行操作，當(dāng)然這也很符合原文的操作。

UMI
UMI就是Unique Molecular Identifier，由4-10個隨機核苷酸組成，在mRNA反轉(zhuǎn)錄后，進入到文庫中，每一個mRNA隨機連上一個UMI，根據(jù)PCR結(jié)果可以計數(shù)不同的UMI，最終統(tǒng)計mRNA的數(shù)量。它的主要作用是，處理PCR 擴增偏差，因為起始文庫很小時需要的PCR擴增次數(shù)就越多，因為越容易引入擴增誤差。

fastq文件更名
為什么要更名？CellRanger 定量過程輸入文件指定命名格式。
如何更名？下圖格式：

image
# 比如，將原來的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 流程

它主要包括四個主要基因表達分析流程：

  • cellranger mkfastq ： 它借鑒了Illumina的bcl2fastq ，可以將一個或多個lane中的混樣測序樣本按照index標(biāo)簽生成樣本對應(yīng)的fastq文件。即對下機數(shù)據(jù)base calling files轉(zhuǎn)為fastq文件。
  • cellranger count ：利用mkfastq生成的fq文件，進行比對(基于STAR)、過濾、UMI計數(shù)。利用細胞的barcode生成gene-barcode矩陣，然后進行樣本分群、基因表達分析
  • cellranger aggr ：接受cellranger count的輸出數(shù)據(jù)，將同一組的不同測序樣本的表達矩陣整合在一起，比如tumor組原來有4個樣本，PBMC組有兩個樣本，現(xiàn)在可以使用aggr生成最后的tumor和PBMC兩個矩陣，并且進行標(biāo)準(zhǔn)化去掉測序深度的影響
  • cellranger reanalyze ：接受cellranger count或cellranger aggr生成的gene-barcode矩陣，使用不同的參數(shù)進行降維、聚類。屬于定制化分析。

它的結(jié)果主要是包含有細胞信息的BAM, MEX, CSV, HDF5 and HTML文件

CellRanger 軟件安裝及測試實戰(zhàn)操作參考：CellRanger走起(四) CellRanger流程概覽。包括安裝詳細指南及其他一些小技巧。

拆分數(shù)據(jù) mkfastq、細胞定量 count、定量組合 aggr

mkfastq 拆分
目的：將每個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指的是下機的原始BCL文件目錄
# 重要的就是測序lane、樣本名稱、index等信息

參考文章CellRanger走起(四) CellRanger流程概覽 中解釋了幾種不同方式輸出fastq文件的不同存放目錄，可根據(jù)實際分析自行查找。

count 定量
這個過程是最重要的，它完成細胞與基因的定量，它將比對、質(zhì)控、定量都包裝了起來，內(nèi)部流程很多，但使用很簡單。每個版本要求的參數(shù)是不同的，尤其是V2與V3版本存在較大差異，這里先對V2進行了解。

# 這是示例，不是真實數(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)的細胞數(shù)量，這個要和實驗設(shè)計結(jié)合起來
# nosecondary 只獲得表達矩陣，不進行后續(xù)的降維、聚類和可視化分析(因為后期會自行用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：比對文件
  • 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：下面進行aggregate使用的文件
  • cloupe.cloupe：官方可視化工具Loupe Cell Browser 輸入文件

count 定量的重點和難點在于參考基因組索引的構(gòu)建，可詳細參考CellRanger走起(四) CellRanger流程概覽，文章中還介紹了怎樣自行構(gòu)建參考基因組索引，及一比對過程的重要細節(jié)信息。

aggr
當(dāng)處理多個生物學(xué)樣本或者一個樣本存在多個重復(fù)/文庫時，最好的操作就是先分別對每個文庫進行單獨的count定量，然后將定量結(jié)果利用aggr組合起來.

  1. 得到count結(jié)果
  2. 構(gòu)建Aggregation CSV
  3. 運行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這個目錄中

cellranger count的結(jié)果

count結(jié)果一般放在out目錄下，主要有summary和analysis兩大類

Summary
是一個HTML格式文件，包括實驗捕獲的細胞數(shù)目、檢測到的基因數(shù)目、測序數(shù)據(jù)的產(chǎn)出與質(zhì)量統(tǒng)計、參考基因組的比對情況。

image

幾個指標(biāo)可以關(guān)注一下：

  • 左上部分中，包括了reads數(shù)、barcodes數(shù)、UMI、index、Q30等統(tǒng)計值

  • 左下是reads比對的比例，包括基因間區(qū)、外顯子、內(nèi)含子，如果比對率太低(一般認為外顯子的比對率要在60%以上)

  • 右上圖是利用barcodes上的UMI標(biāo)簽分布來估計細胞數(shù)，綠/藍色表示細胞，灰色表示背景，其中Y軸表示每個barcode對應(yīng)的UMI數(shù)量，X軸是一定數(shù)量的UMI序列所對應(yīng)barcode的數(shù)量，比如上圖中有1000個barcodes含有10k個UMI，細胞過濾就是通過這個圖來展示的。

  • 首先明確，barcode用來區(qū)分細胞，UMI用來區(qū)分轉(zhuǎn)錄本；其次，barcodes數(shù)量時要大于細胞數(shù)量的（以保證每個細胞都會有barcode來進行區(qū)分）如果圖線陡降說明

另參考文章CellRanger走起(二) CellRanger out 結(jié)果 中還記錄了一個手動構(gòu)建物種GTF文件的案例嘗試，及相應(yīng)的檢驗方法，可以根據(jù)需要自行學(xué)習(xí)。

總之，無論是從公共數(shù)據(jù)庫下載單細胞SRA數(shù)據(jù)，還是來自10x 平臺的原始下機數(shù)據(jù)，或者先經(jīng)過fastq-dump轉(zhuǎn)為fastq文件，然后質(zhì)控，然后CellRanger count 定量，或者首先自行拆分數(shù)據(jù)mkfastq，然后構(gòu)建索引進行細胞定量，及可選的定量組合等。CellRanger的總體流程見本文記錄，但是其中的各個細節(jié)仍需要自行實踐總結(jié)。

作者：新欣enjoy
鏈接：http://www.itdecent.cn/p/5157ab9f6977

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容