scRNA-seq表達矩陣的構建

• 2020 年 3 月 27 日
• 筆記

希望大家能有所收穫！

目錄

⊙引言—關於課程

⊙scRNA-seq簡介

⊙scRNA-seq原始數據的質控

⊙scRNA-seq數據處理—文件格式小結

⊙scRNA-seq數據處理—demultiplexing

⊙scRNA-seq數據的處理—STAR

⊙scRNA-seq數據處理—Kallisto

正文

表達矩陣的構建

scRNA-seq數據的許多分析以表達矩陣為起點。按照慣例，表達矩陣的每一行代表一個基因，每列代表一個細胞（儘管一些作者使用轉置矩陣）。每個條目代表給定細胞中特定基因的表達水平。基因表達的測量單位取決於protocol和使用的一般方式。

4.1 reads QC

scRNA-seq實驗的結果是大量的cDNA reads的集合。第一步是確保reads具有高品質。可以使用一些常規工具（如FastQC或Kraken）執行品質控制。

假設我們的reads位於experiment.bam中，我們將運行FastQC

$<path_to_fastQC>/fastQC experiment.bam

下面是FastQC輸出125 bp reads數據集的示例。該圖顯示了一個技術錯誤，導致幾個鹼基無法在read中心無法正確讀取。但是，由於讀取的其餘部分品質很高，因此該錯誤很可能對比對效率的影響可以忽略不計。

此外，使用Integrative Genomics Browser（IGV）或SeqMonk對數據進行可視化通常很有幫助。

4.2 reads比對

在從讀數中trim掉低品質鹼基後，可以將剩餘序列比對到參考基因組。同樣，不需要特殊的方法，因此我們可以使用STAR或TopHat比對工具。對於有良好注釋基因組的生物（例如小鼠，人）的全轉錄數據集，偽比對方法（例如Kallisto，Salmon）可能優於常規比對。對於具有數十或數十萬個reads的基於drop-seq的數據集，偽比對可能更有效，因為它們的運行時間可能比傳統的比對工具低幾個數量級。

使用STAR來比對reads.bam的例子：

$<path_to_STAR>/STAR --runThreadN 1 --runMode alignReads  --readFilesIn reads1.fq.gz reads2.fq.gz --readFilesCommand zcat --genomeDir <path>  --parametersFiles FileOfMoreParameters.txt --outFileNamePrefix <outpath>/output

注意，如果使用spike-ins，參考序列應該與所述的DNA序列spike-ins在比對之前的分子。

請注意，使用UMI時，應從讀取序列中刪除其barcode。通常的做法是將barcode添加到read名稱上。

一旦將每個細胞的reads比對到參考基因組，我們需要確保每個細胞的足夠數量的reads可以比對到參考基因組。根據我們的經驗，小鼠或人類細胞的可比對的reads比例為60-70％。但是，此結果可能會因protocol，reads長度和reads比對的設置而異。一般來說，我們希望所有細胞都具有相似的比對的reads部分，因此應檢查並可能刪除任何異常值。低比例的可比對reads通常意味著污染。

如何使用Salmon量化表達的一個例子是

$<path_to_Salmon>/salmon quant -i salmon_transcript_index -1 reads1.fq.gz -2 reads2.fq.gz -p #threads -l A -g genome.gtf --seqBias --gcBias --posBias

注意 Salmon根據我們的經驗產生估計reads數和估計的每百萬轉錄數（tpm），後者用於校正scRNASeq的長基因的表達，因此我們建議使用reads數。

4.3 比對例子

下面的直方圖顯示了scRNA-seq實驗中映射到每個細胞的讀數總數。每個條形代表一個細胞，它們按照每個細胞的總讀數按升序排序。三個紅色箭頭表示在覆蓋範圍方面為異常值的細胞，應將其從進一步分析中刪除。兩個黃色箭頭指向具有令人驚訝的大量未映射讀數的細胞。在該實施例中，我們在比對QC步驟期間保持細胞，但是由於核糖體RNA讀取的高比例，它們隨後在細胞QC期間被去除

4.4 對比QC

在將原始測序映射到基因組後，我們需要評估映射的品質。有許多方法可以測量繪圖品質，包括：映射到rRNA / tRNA的讀數量，獨特映射讀數的比例，跨剪接點的讀數映射，沿著轉錄本的讀取深度。為大量RNA-seq開發的方法，如RSeQC，適用於單細胞數據：

python <RSeQCpath>/geneBody_coverage.py -i input.bam -r genome.bed -o output.txt  python <RSeQCpath>/bam_stat.py -i input.bam -r genome.bed -o output.txt  python <RSeQCpath>/split_bam.py -i input.bam -r rRNAmask.bed -o output.txt

然而，預期結果將取決於實驗方案，例如許多scRNA-seq方法使用poly-A選擇以避免對rRNA進行測序，這導致跨基因的讀取覆蓋中的3'偏差（也稱為基因體覆蓋）。下圖顯示了這個3'偏差以及三個單元格，它們是異常值並從數據集中刪除：

4.5 reads量化

下一步是量化每個細胞的每個基因的表達水平。對於mRNA數據，我們可以使用為大量RNA-seq數據開發的工具之一，例如HT-seq或FeatureCounts

# include multimapping  <featureCounts_path>/featureCounts -O -M -Q 30 -p -a genome.gtf -o outputfile input.bam  # exclude multimapping  <featureCounts_path>/featureCounts -Q 30 -p -a genome.gtf -o outputfile input.bam

獨特的分子標識符（UMI）使得計算分子的絕對數量成為可能，並且它們已被證明在scRNA-seq中很受歡迎。我們將在下一章討論如何處理UMI。

4.6 獨特的分子標識符（UMI）

感謝EMBL Monterotondo的 Andreas Buness 在本節的合作。

4.6.1 簡介

獨特的分子標記是在反轉錄過程中添加到轉錄本中的短（4-10bp）隨機條形碼。它們使測序讀數能夠分配到單個轉錄物分子，從而從scRNASeq數據中去除擴增雜訊和偏差。

當對包含UMI的數據進行排序時，使用技術僅對包含UMI（通常為3'末端）的轉錄本的末端進行特異性測序。

4.6.2 比對barcode

由於獨特的條形碼數量（4N,其中N是UMI的長度）遠小於每個細胞的分子總數（~106），每個條形碼通常被分配給多個轉錄本。因此，為了識別獨特的分子，必須使用條形碼和映射位置（轉錄物）。第一步是映射UMI讀數，我們建議使用STAR，因為它很快並輸出高品質的BAM比對。此外，映射位置可用於例如。識別記錄不良的3'UTR成績單。

UMI測序通常由配對末端讀數組成，其中每對讀取一個讀取細胞和UMI條形碼，而另一個讀取包含來自轉錄物的外顯子序列（圖4.5）。注意，建議修剪和/或過濾以去除含有poly-A序列的讀段，以避免由於這些讀取映射到具有內部poly-A / poly-T序列的基因/轉錄物而導致的錯誤。

處理UMI實驗中的讀取後，通常使用以下約定：

1. UMI被添加到另一個配對讀取的讀取名稱中。
2. 讀取按單元條形碼分類到單獨的文件中

• 對於極大的淺數據集，可以將單元條形碼添加到讀取名稱中以減少文件數量。

4.6.3 比對barcode

理論上，每個獨特的UMI轉錄物對應代表源自單個RNA分子的所有讀數。但是，實際情況往往並非如此，最常見的原因是：

1. 不同的UMI不一定意味著不同的分子

• 由於PCR或測序錯誤，鹼基對取代事件可導致新的UMI序列。較長的UMI會產生更多出錯的機會，並且根據單元條碼的估計，我們預計7-10％的10bp UMI至少會包含一個錯誤。如果沒有糾正，這種類型的錯誤將導致高估轉錄本的數量。

1. 不同的轉錄物不一定意味著不同的分子

• 映射錯誤和/或多映射讀取可能導致某些UMI被分配給錯誤的基因/轉錄本。這種類型的錯誤也會導致高估轉錄本的數量。

1. 相同的UMI不一定意味著相同的分子

• UMI頻率和短UMI的偏差可導致相同的UMI附著於來自相同基因的不同mRNA分子。因此，可能低估了轉錄本的數量。

4.6.4 糾正誤差

如何最好地解釋UMI中的錯誤仍然是一個活躍的研究領域。我們知道解決上述問題的最佳方法是：

1. UMI工具的定向鄰接方法實現了一個過程，該過程考慮了不匹配的數量和類似UMI的相對頻率，以識別可能的PCR /排序錯誤。
2. 目前是一個未決問題。通過刪除具有少量讀取的UMI來支援它們與特定轉錄本的關聯，或者通過移除所有多映射讀取，可以減輕該問題。
3. 由Grun，Kester和van Oudenaarden（2014）提出的簡單飽和度（又名「碰撞概率」）校正來估計分子的真實數量M：
4. 其中N =唯一UMI條形碼的總數，n =觀察到的條形碼數。

這種方法的一個重要警告是它假設所有UMI都是同等頻繁的。在大多數情況下，這是不正確的，因為通常存在與GC內容相關的偏差。

確定如何最好地處理和使用UMI目前是生物資訊學界的一個活躍的研究領域。我們知道最近開發的幾種方法，包括：

• UMI工具
• PoissonUMIs
• zUMIs
• dropEst

4.6.5 下游分析

當前的UMI平台（DropSeq，InDrop，ICell8）具有低且高度可變的捕獲效率，如下圖所示。

這種可變性會引入強烈的偏差，需要在下游分析中加以考慮。最近的分析通常基於細胞類型或生物途徑將細胞/基因聚集在一起以增加功率。對這些數據的穩健統計分析仍然是一個開放的研究問題，還有待確定如何最好地調整偏差。

練習1我們為您提供了由三個不同個體產生的誘導多能幹細胞的UMI計數和讀數（Tung et al.2017）（有關此數據集的詳細資訊，請參閱：第7.1章）。

umi_counts <- read.table("tung/molecules.txt", sep = "t")  read_counts <- read.table("tung/reads.txt", sep = "t")

使用此數據：

1. 繪製捕獲效率的可變性
2. 確定擴增率：每個UMI的平均讀數。