CellRanger原理

0. 概述

Cell Ranger 是 10X genomics 官方提供的一套針對單細胞 RNA 測序輸出結(jié)果進行比對、定量、聚類及基因表達分析的分析流程,它包含有與單細胞基因表達分析相關(guān)的四個pipelines,分別是:
cellranger mkfastq 流程:其功能為將 Illumina 測序儀產(chǎn)生的 raw base call (BCL) 文件解析成 FASTQ 文件。

cellranger count 流程:其功能為將 cellranger mkfastq 產(chǎn)生的或其他來源的 FASTQ 文件進行比對、過濾、barcode 計數(shù)以及 UMI 計數(shù),并可以生成 feature-barcode 定量矩陣,隨后確定細胞群并進行基因表達分析。

cellranger aggr 流程:接受cellranger count的輸出數(shù)據(jù),將同一組的不同測序樣本的表達矩陣整合在一起,比如tumor組原來有4個樣本,PBMC組有兩個樣本,現(xiàn)在可以使用aggr生成最后的tumor和PBMC兩個矩陣,并且進行標準化去掉測序深度的影響。

cellranger reanalyze 流程:接受 cellranger count 或 cellranger aggr 的輸出文件,對數(shù)據(jù)重新進行降維、聚類等后續(xù)分析。

以上四個pipeline 均將轉(zhuǎn)錄組常用比對軟件 STAR 封裝其中,可以輸出帶有細胞信息的 BAM、MEX、CSV、HDF5 及 HTML 等格式的結(jié)果。

參考:CellRanger官網(wǎng)

1. mkfastq 流程

如果是直接拿到的R1/R2的fastq文件,那么就直接上cellranger count。

illumina下機fastq文件命名規(guī)則:

但是如果是I1/R1/R2的數(shù)據(jù),那就還要跑個cellranger mkfastq。

2. count流程

2.1命令
/software/cellranger-4.0.0/cellranger count \
--id=cellranger_HH \ #數(shù)據(jù)結(jié)果文件夾的名稱
--fastqs=rawdata \ #原始文件路徑
--sample=Con-1-10X5 \ # 原始下機數(shù)據(jù)文件的前綴 (fastq.gz文件名當(dāng)中S1之前的字段)
--transcriptome=/reference/refdata-gex-GRCh38-2020-A #參考基因組文件目錄路徑
--description=H2022-xxxx_人-PBMC \ # 對樣本信息進行描述
--localcores=16 \ # 將限制 cellranger 一次最多使用 16 個核
--localmem=64 \ # 將限制 cellranger 使用的最大內(nèi)存量 64(GB)

注意:
> --localcores用于設(shè)置程序的核心數(shù),不設(shè)置的話,默認是使用系統(tǒng)所能夠使用的最多的核心。(設(shè)置是為了避免cellranger耗盡所有資源)
> --localmem :限制cellrange使用的內(nèi)存數(shù)(設(shè)置是為了避免cellranger耗盡所有資源)
> 如果是裂核項目,需加上--include-introns 參數(shù)運行。因為添加此參數(shù)可以將 reads 映射到內(nèi)含子區(qū)域,提高對含有大量前 mRNA 分子(如細胞核)樣品的敏感性。
> 如果遇到結(jié)果不符合預(yù)期,可根據(jù)情況利用 --force-cells=Num 指定細胞數(shù),通過調(diào)整細胞數(shù)以獲得預(yù)期的結(jié)果。
> --except-cells:期望得到的細胞數(shù)目,默認是3000個,一般大家都設(shè)置1000

最好不要去使用--except-cells--force-cells=Num這兩個參數(shù),沒必要去人為參與細胞識別的過程。

2.2 結(jié)果文件

我們主要關(guān)注outs這個文件夾。outs文件夾中的結(jié)果文件主要如下

Output Files:
Outs
├── analysis # 數(shù)據(jù)分析結(jié)果文件
│   ├── clustering # 聚類文件夾,圖聚類,k-means 聚類
│   ├── diffexp # 差異分析
│   ├── pca # pca 線性降維
│   ├── tsne # tsne 非線性降維
│   └── umap # umap 非線性降維
├── cloupe.cloupe # 用于 Loupe Browser 可視化和分析的輸入文件
├── filtered_feature_bc_matrix # 過濾后的 barcode 矩陣信息,用于后續(xù)分析的必要文件│   
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5 # 過濾后的 barcode 信息 HDF5 格式
├── metrics_summary.csv # 結(jié)果匯總 csv 文件
├── molecule_info.h5 # 實驗相關(guān)的文庫,GEM,barcode表達量,UMI 等信息的HDF5格式的文件,cellranger aggr aggregate 數(shù)據(jù)的時候會用到的文件
├── possorted_genome_bam.bam # 比對到參考基因組和轉(zhuǎn)錄本上并帶有 barcode 信息的 reads 信息
├── possorted_genome_bam.bam.bai # possorted_genome_bam.bam 的索引文件
├── raw_feature_bc_matrix # 未經(jīng)過濾的所有 barcode 的矩陣信息
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5 # 原始 barcode 信息 HDF5 格式
└── web_summary.html # 結(jié)果匯總 html 文件及可視化

3. out結(jié)果文件解讀

3.1Unfiltered feature-barcode matricesFiltered feature-barcode matrices是怎么來的

這兩個文件中都包括了如下三個文件,Unfiltered目錄下是所有的barcode信息,包含了細胞相關(guān)的barcoed和背景barcode,而Filtered目錄下只包含細胞相關(guān)的barcode信息。下游分析我們需要的就是filtered_gene_bc_matrices 文件夾下的文件作為輸入文件輸入到seurat或者scanpy。

├── barcodes.tsv
├── genes.tsv
└── matrix.mtx
  • Unfiltered feature-barcode matrices

(1)修剪Reads:針對 3’ 建庫數(shù)據(jù)的基因表達比對,在比對之前會先對 reads 進行修剪。
cDNA 的全長結(jié)構(gòu)中,在 3’ 和 5’ 端分別帶有 poly-A 尾和TSO 序列結(jié)構(gòu)(相對于比較長的 RNA 分子,一部分來自短 RNA 分子的 reads 可能僅包含 TSO 和 poly-A 序列的其中一種)。由于這種低復(fù)雜度的非模板序列的存在有可能混淆 reads 的映射,所以在比對之前一般會將 poly-A 尾和 TSO 序列分別從 reads 的 3’ 端和 5’ 端切除,這一步驟有助于提高分析的靈敏度和軟件分析的效率。
(2)判斷 reads 是否比對到了基因組
Cell Ranger 中封裝了比對軟件 STAR,根據(jù)轉(zhuǎn)錄本的注釋文件 GTF 中的注釋信息,使用 STAR 來判斷reads 是比對到了外顯子、內(nèi)含子還是基因間區(qū)上,或者說來判斷 reads 是否比對到了基因組上。
當(dāng)一條 read 至少要有 50% 堿基序列與基因組上的外顯子堿基互補配對,認為其比對到了外顯子上;若 reads 未比對上外顯子但與內(nèi)含子相交,則認為其比對到了內(nèi)含子上;否則為比對到了基因間區(qū)。若 reads 比對到了一個單一的外顯子位點,但同時比對到了一個或多個非外顯子位點,則優(yōu)先認為該 read 比對到了外顯子位點,MAPQ 為 255。
(3)判斷 reads 是否比對到了轉(zhuǎn)錄本
Cell Ranger 通過檢測 reads 比對上的外顯子和內(nèi)含子與轉(zhuǎn)錄本的相容性,進一步將 reads 與注釋的轉(zhuǎn)錄本對齊。如下圖所示,reads 根據(jù)它們是正義還是反義,以及它們是外顯子還是內(nèi)含子,或者它們的剪接模式是否與該基因相關(guān)的轉(zhuǎn)錄本注釋兼容來分類。

上圖中,綠色展示的是基因及基因中所包含的外顯子,Transcript 1 和 Transcript 2 為基因經(jīng)過可變剪切形成的兩種轉(zhuǎn)錄本所包含的外顯子。針對比對到正義鏈上的reads,如果 reads 比對到了一個外顯子上或者比對到兩個相鄰的外顯子上,則該 read 被分類為轉(zhuǎn)錄本 read(藍色);如果 reads 比對到兩個不相鄰的外顯子上,則該 read 被分類為外顯子 read(淺藍色);如果 reads 比對到內(nèi)含子區(qū)域,則該 read 被分類為內(nèi)含子 read(紅色);紫色表示 reads 比對到反義鏈上。

小知識(單細胞測序和單核細胞測序的區(qū)別)
在默認情況下,只有藍色的轉(zhuǎn)錄本 read 會被計入到 UMI 計數(shù)中。但在某些情況下,如在實驗時輸入的為細胞核時,未剪接的轉(zhuǎn)錄本有可能產(chǎn)生高水平的內(nèi)含子序列,為了將這些內(nèi)含子 read 計入,cellranger count 可以添加一個參數(shù)為 include-introns。當(dāng)使用該參數(shù)時,任何比對到單個基因的 reads ---- 包括轉(zhuǎn)錄本 read(藍色)、外顯子 read(淺藍色)和內(nèi)含子 read(紅色)都會計入 UMI 計數(shù)中。
此外,只有在基因組上有唯一比對位點的 reads 才被計入到UMI計數(shù)中。

(4)進行 UMI 計數(shù)
在計算 UMIs 之前,Cell Ranger 會試圖矯正 UMI 序列中的測序錯誤。
在轉(zhuǎn)錄本上有唯一比對位點的 reads 根據(jù)他們的barcode、UMI 和比對到的基因被分成不同的組。如果兩個組的 reads 擁有相同的 barcode 序列并比對到同一個基因上,但是 UMI 序列中有一個堿基不同,那么其中一個 UMI 有可能是因為測序中的堿基替換錯誤而引入的。在這種情況下,UMI 的reads 數(shù)量少的那一組會被更正為 UMI 的reads數(shù)量多的那組。
矯正可能的測序錯誤后進行 UMI 計數(shù)。
Cell Ranger 會再次按照 UMI(可能是修正后的)、barcode 和比對到的基因?qū)?reads 進行分組。如果兩組或者多組的 reads 擁有相同的 barcods 和 UMI 序列,但是比對到了不同的基因上,那么 reads 計數(shù)最高的那組比對到的基因會被進行 UMI 計數(shù),其他的組則被舍棄掉。如果 reads 最高計數(shù)相同,則全部的組都被舍棄掉。

經(jīng)過這兩步過濾步驟后,每一個被統(tǒng)計到的barcode、UMI 和 基因都會被保存在未過濾的 feature-barcode 矩陣中,輸出在 unfiltered feature-barcode matrix 文件夾中。

  • Filtered feature-barcode matrices

Cell Ranger 進行細胞識別的算法直接影響了我們可以獲得的高質(zhì)量細胞的數(shù)量及數(shù)據(jù)質(zhì)量。在基于 droplet 方法的單細胞技術(shù)中,通常認為含有細胞的液滴應(yīng)該含有更多的RNA,因此其在 UMI 總量上應(yīng)該與空細胞(背景噪音)存在明顯的區(qū)分(也就是我們常說的 barcodes 排序圖上的拐點)。然而實際上微滴間不同的擴增效率會導(dǎo)致一些較小的細胞跟空細胞在 UMI 總數(shù)上是相近的,無法僅通過 UMI 總數(shù)很好地區(qū)分空細胞和非空細胞。尤其當(dāng)樣本中混雜了不同大小的細胞,例如在腫瘤樣本中一般混雜著體型較大的腫瘤細胞和較小的腫瘤浸潤淋巴細胞,浸潤淋巴細胞則較難與空細胞區(qū)分。

為解決這一難題,Cell Ranger 的算法采用了兩步法分別基于 UMI 閾值識別高 RNA 含量細胞以及基于表達譜識別低 RNA 含量細胞
(1)第一步,選取一個 UMI 總數(shù)的閾值,所有大于這個閾值的 barcodes 被識別為細胞。這一步保證了高 RNA 含量的 barcodes 被保留。
具體的算法是:將 UMI 計數(shù)從高到低進行排序,根據(jù)預(yù)期細胞數(shù) N(軟件默認N=3000),排名前 N 個細胞中的 99 分位 UMI 數(shù)值記為 m,將所有 UMI 計數(shù)大于 m/10 這一閾值的 barcodes 標記為高質(zhì)量細胞。
(2)剩余未通過閾值的 barcodes 則進行第二步的篩選,根據(jù)與空細胞 RNA 表達譜是否存在顯著差異來回收潛在的低 RNA 含量細胞。
此算法基于 Lun et al. 2019 年發(fā)表的算法 EmptyDrops。

a. 首先選取一組低 UMI 計數(shù)的 barcodes 作為背景集(來代表空細胞),用這部分 barcodes 的表達譜構(gòu)建一個“背景模型”。

b. 接著將在第一步驟中所有未被標注為高質(zhì)量細胞的 barcodes 依次和背景模型相比較,那些與背景模型存在顯著差異的細胞會被回收到高質(zhì)量細胞的行列。

例如:
下圖是一群高 RNA 含量的 293T細胞和一群低 RNA 含量的 PBMC 細胞的混合樣本數(shù)據(jù)??梢钥吹皆诒粯擞洖楦哔|(zhì)量的部分出現(xiàn)了兩個群體(由第一個拐點A 大致分開),在第二個拐點 B 附近的區(qū)域同時包含空細胞和高質(zhì)量細胞,這部分細胞即為從第二步中回收出的細胞,圖片中顏色的深淺代表了局部高質(zhì)量細胞的比例。

所有被回收的高質(zhì)量細胞的矩陣,會被輸出到 filtered_feature_bc_matrix 文件夾中,根據(jù)矩陣信息進行下游分析。

3.2 網(wǎng)頁報告解讀

網(wǎng)頁的結(jié)果分成了summary和analysis兩部分, summary部分包含如下結(jié)果:

參考:
Cell Ranger 知多少?(上)
Cell Ranger 知多少?(中)
Cell Ranger 知多少(下)
【單細胞測序】 關(guān)于cellranger

最后編輯于
?著作權(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ù)。
禁止轉(zhuǎn)載,如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

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

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