Hemberg-lab單細胞轉錄組數據分析（一）

Hemberg-lab單細胞轉錄組數據分析（二）

Hemberg-lab單細胞轉錄組數據分析（三）

Hemberg-lab單細胞轉錄組數據分析（四）

Hemberg-lab單細胞轉錄組數據分析（五）

收藏|北大生信平臺"單細胞分析、染色質分析"視頻和PPT分享

該如何自學入門生物信息學

生物信息之程序學習

構建表達矩陣

scRNA-seq數據的許多分析以表達矩陣為起點。一般來講，表達矩陣的每一行代表一個基因，每一列代表一個細胞（但是一些作者會做個轉置）。每個條目代表特定基因在給定細胞中的表達水平。而表達值的測量單位取決于建庫方案和所用的標準化方法。

reads質控

見前面章節FastQC部分。

另外，使用Integrative Genomics Browser（IGV）或SeqMonk通常對數據可視化很有幫助，具體見下。

• 測序數據可視化 (一)

• IGV基因組瀏覽器可視化高通量測序數據

• 高通量數據分析必備-基因組瀏覽器使用介紹 - 1

• 高通量數據分析必備-基因組瀏覽器使用介紹 - 2

• 高通量數據分析必備-基因組瀏覽器使用介紹 - 3

Reads比對

見前面章節STAR部分和Kallisto部分。

注釋的很好的模式生物（例如小鼠，人）有著大量全長轉錄本數據集，偽比對方法（例如Kallisto，Salmon)可能優于常規比對方法。drop-seq方法獲得的數據集有數以千萬條reads，偽比對工具的運行時間比傳統比對工具會少幾個數量級，更有時間優勢。從39個轉錄組分析工具，120種組合評估(轉錄組分析工具哪家強-導讀版)一文中可以看出，偽比對工具的準確性和穩定性也相對比較高。

用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（已知濃度的外源RNA分子），在比對前應該將參考基因組和spike-in分子的DNA序列合并作為共同”參考基因組”。具體見之前文件格式部分。

注意，使用UMI時，應從read序列中刪除其條形碼。常見的是將條形碼加到read名稱上。

一旦reads完成了到基因組的比對，我們需要檢查比對率和確保有足夠多的reads比對回了參考基因組。根據我們的經驗，小鼠或人類細胞中read的比對率為60-70％。但是這個結果可能會因建庫方案、read長度和比對工具參數設置而有所不同。常規上，我們希望所有細胞都具有相似的read比對率。如果有樣品比對率異常低或比對回去的reads異常低，則需要多加注意甚至從后續分析中移除。較低的read比對率通常表示存在污染。生信寶典建議取出幾十條未必對回去的reads做個blast，看下能否比對到其它物種來確定是污染還是測序錯誤還是比對參數設置的問題。

一個用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操作會得到一個估計的read counts和transcripts per million(tpm)。根據我們的經驗，TPM對單細胞測序中長基因的表達做了過度校正，因此我們建議使用read counts。

比對示例

下面直方圖 (http://www.ehbio.com/ImageGP)顯示了scRNA-seq實驗中每個細胞不同比對狀態的reads的數目。每個柱子代表一個細胞，按細胞的總read數升序排列。三個紅色箭頭標記的是比對到基因組的reads較低的異常樣本，應該在后續分析中移除。兩個黃色箭頭指的是unmapped reads數目十分大的細胞。該例中，在比對質控期間這兩個細胞會保留下來，但后期細胞質控時這兩個細胞會因為核糖體RNA?reads比例過高而移除。

Mapping QC

在把原始序列比對到基因組后，需要評估比對質量。這可以從多個角度進行評估，包括：rRNA/tRNAs的reads的占比或總量，reads在基因組上唯一比對位置的比例，比對到splice junction的reads比例，reads在轉錄本的覆蓋均一性或深度。而為Bulk RNA-seq開發的方法如RSeQC，也適用于單細胞數據：

#pip install RSeQC geneBody_coverage.py -i input.bam -r genome.bed -o output.txt bam_stat.py -i input.bam -r genome.bed -o output.txt split_bam.py -i input.bam -r rRNAmask.bed -o output.txt

然而，預期結果的評估取決于采用的建庫方案，例如許多scRNA-seq用poly-A selection捕獲轉錄本。這個方法可以排除核糖體RNA污染，但會導致3'區域更容易測到。下圖展示了測序reads分布的3'偏好性，和去除的三個異常細胞的結果 (應該是最下面3條，推測是降解嚴重)。

Reads量化

scRNA-seq基因定量計算可以用bulk RNA-seq一樣的工具，比如HT-seq or FeatureCounts。

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

唯一分子標識符UMI讓計算轉錄本的絕對量成為可能，并且在scRNA-seq中很受歡迎。我們將在下一章討論如何處理UMI。

唯一分子標識符（Unique Molecular Identifiers, UMI）

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

UMI添加到每個轉錄本上

唯一分子標識符 (UMI)是在反轉錄過程中添加到轉錄本上的短的（4-10 bp）隨機條形碼序列。它們使得測序reads可以對應到單個轉錄本，從而去除擴增噪聲和偏好性。

當測序含UMI的文庫時，僅對包含UMI的轉錄本末端 (通常為3'末端)進行性測序。

比對UMI條形碼

由于UMI數量（, N是UMIs的長度值）比每個細胞中的RNA分子數（~）少得多，每個UMI條形碼可能會連接到多個轉錄本，因此需要借助條形碼序列和reads比對位置兩個條件鑒定起始的轉錄本分子。第一步是比對UMI reads，推薦用STAR來處理，因為它處理速度快且輸出高質量的BAM比對。此外，比對位置的準確性對識別新的3'UTR區域很有意義。

UMI測序通常由雙端reads組成，其中一端read是捕獲細胞和UMI的條形碼，而另一端read包含轉錄本的外顯子序列。注意，推薦去除reads中的poly-A序列部分，以免這些reads比對到轉錄本內部poly-A或poly-T序列而產生錯誤。

處理UMI實驗中的reads，通常有以下慣例：

UMI被添加到另一個配對read的序列名稱中。

reads按細胞條形碼分類到單獨的文件中 (見前面的文章)。但對于細胞量極大的低深度測序數據集 (drop-seq)，可以將細胞條形碼添加到read名稱中而不是拆分為單獨文件以減少文件數量。

Counting 條形碼

理論上，每個唯一的UMI-轉錄本對應該對應來源于一個RNA分子的所有reads。但是現實往往并非如此，最常見的原因是：

不同的UMI序列不一定表示它們是不同的UMI分子由于PCR或測序錯誤，堿基替換可能導致新的UMI序列。較長的UMI出現堿基替換的機會更多。根據細胞條碼測序錯誤估計，7-10％的10 bp長度的UMI中至少有一個堿基替換。如果錯誤沒有糾正，將會過高估計轉錄本的數目。

不同的轉錄本不一定是不同的分子比對錯誤或多個比對位置可能導致某些UMI對應到錯誤的基因/轉錄本。這種類型的錯誤也會導致過高估計轉錄本的數目。

相同的UMI不一定意味著相同的分子UMI頻次的不同和短UMI可導致同一UMI和相同基因的不同mRNA分子相連，進而可能低估轉錄本數量。

錯誤糾正

如何最好的校正UMIs中的這些問題仍然是一個活躍的研究領域。我們自己認為的最好的解決上述問題的方法是：

UMI-tools’，設計了directional-adjacency算法，同時考慮錯配數目和相似UMI的相對頻率來識別PCR和測序錯誤。(alevin, cellranger都是不錯的選擇，后面詳細介紹)

問題還無法完全解決，但通過刪除只有很少read支持的UMIs-轉錄本對，或者移除所有多比對位置的reads，可能會減輕該問題。

Simple saturation (也稱為”collision probability”)方法來估計分子的數量

其中N=唯一UMI條形碼的總數，n=觀察到的條形碼數。

這個方法的一個重要缺陷是它假設所有UMI出現頻率相同。但因為序列GC含量不同引入的偏差使得這一假設在大多數情況下這是不正確的。

如何最好地處理和使用UMI在目前生物信息學界是一個活躍的研究領域。而我們了解到的幾種最近開發的方法有：

• UMI-tools

• PoissonUMIs

• zUMIs

• dropEst

下游分析

當前的UMI平臺（DropSeq，InDrop，ICell8）展現出從低到高變化很大的捕獲效率，如下圖所示。

這一高可變性可能會引入很強的偏差，需要在下游分析時考慮到。現在的分析通常根據細胞類型或生物通路把細胞/gene混合一起增加檢測能力。更合適的統計分析方法亟待研究以便更好地調整這些偏差，使得結果更能反映真實現象。

練習1?數據是三個不同來源的誘導多功能干細胞的UMI counts和read counts?(有關此數據集的詳細信息請參閱后續文章)。

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

使用此數據：

繪制捕獲效率的變化

確定擴增率：每個UMI對應的平均reads數。

# Exercise 1 # Part 1 plot(colSums(umi_counts), colSums(umi_counts > 0), xlab="Total Molecules Detected", ylab="Total Genes Detected")# Part 2 amp_rate <- sum(read_counts)/sum(umi_counts) amp_rate

總結

以上是生活随笔為你收集整理的Hemberg-lab单细胞转录组数据分析（六）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。