Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（一）

Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（二）

Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（三）

Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（四）

Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（五）

收藏|北大生信平臺(tái)"單細(xì)胞分析、染色質(zhì)分析"視頻和PPT分享

該如何自學(xué)入門生物信息學(xué)

生物信息之程序?qū)W習(xí)

構(gòu)建表達(dá)矩陣

scRNA-seq數(shù)據(jù)的許多分析以表達(dá)矩陣為起點(diǎn)。一般來講，表達(dá)矩陣的每一行代表一個(gè)基因，每一列代表一個(gè)細(xì)胞（但是一些作者會(huì)做個(gè)轉(zhuǎn)置）。每個(gè)條目代表特定基因在給定細(xì)胞中的表達(dá)水平。而表達(dá)值的測量單位取決于建庫方案和所用的標(biāo)準(zhǔn)化方法。

reads質(zhì)控

見前面章節(jié)FastQC部分。

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

• 測序數(shù)據(jù)可視化 (一)

• IGV基因組瀏覽器可視化高通量測序數(shù)據(jù)

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

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

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

Reads比對(duì)

見前面章節(jié)STAR部分和Kallisto部分。

注釋的很好的模式生物（例如小鼠，人）有著大量全長轉(zhuǎn)錄本數(shù)據(jù)集，偽比對(duì)方法（例如Kallisto，Salmon)可能優(yōu)于常規(guī)比對(duì)方法。drop-seq方法獲得的數(shù)據(jù)集有數(shù)以千萬條reads，偽比對(duì)工具的運(yùn)行時(shí)間比傳統(tǒng)比對(duì)工具會(huì)少幾個(gè)數(shù)量級(jí)，更有時(shí)間優(yōu)勢。從39個(gè)轉(zhuǎn)錄組分析工具，120種組合評(píng)估(轉(zhuǎn)錄組分析工具哪家強(qiáng)-導(dǎo)讀版)一文中可以看出，偽比對(duì)工具的準(zhǔn)確性和穩(wěn)定性也相對(duì)比較高。

用STAR比對(duì)的操作示例 (前面章節(jié)部分更詳細(xì))

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分子），在比對(duì)前應(yīng)該將參考基因組和spike-in分子的DNA序列合并作為共同”參考基因組”。具體見之前文件格式部分。

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

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

一個(gè)用Salmon量化表達(dá)操作的例子

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操作會(huì)得到一個(gè)估計(jì)的read counts和transcripts per million(tpm)。根據(jù)我們的經(jīng)驗(yàn)，TPM對(duì)單細(xì)胞測序中長基因的表達(dá)做了過度校正，因此我們建議使用read counts。

比對(duì)示例

下面直方圖 (http://www.ehbio.com/ImageGP)顯示了scRNA-seq實(shí)驗(yàn)中每個(gè)細(xì)胞不同比對(duì)狀態(tài)的reads的數(shù)目。每個(gè)柱子代表一個(gè)細(xì)胞，按細(xì)胞的總read數(shù)升序排列。三個(gè)紅色箭頭標(biāo)記的是比對(duì)到基因組的reads較低的異常樣本，應(yīng)該在后續(xù)分析中移除。兩個(gè)黃色箭頭指的是unmapped reads數(shù)目十分大的細(xì)胞。該例中，在比對(duì)質(zhì)控期間這兩個(gè)細(xì)胞會(huì)保留下來，但后期細(xì)胞質(zhì)控時(shí)這兩個(gè)細(xì)胞會(huì)因?yàn)楹颂求wRNA?reads比例過高而移除。

Mapping QC

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

#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

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

Reads量化

scRNA-seq基因定量計(jì)算可以用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

唯一分子標(biāo)識(shí)符UMI讓計(jì)算轉(zhuǎn)錄本的絕對(duì)量成為可能，并且在scRNA-seq中很受歡迎。我們將在下一章討論如何處理UMI。

唯一分子標(biāo)識(shí)符（Unique Molecular Identifiers, UMI）

感謝EMBL Monterotondo的 Andreas Buness 在本節(jié)的合作。

UMI添加到每個(gè)轉(zhuǎn)錄本上

唯一分子標(biāo)識(shí)符 (UMI)是在反轉(zhuǎn)錄過程中添加到轉(zhuǎn)錄本上的短的（4-10 bp）隨機(jī)條形碼序列。它們使得測序reads可以對(duì)應(yīng)到單個(gè)轉(zhuǎn)錄本，從而去除擴(kuò)增噪聲和偏好性。

當(dāng)測序含UMI的文庫時(shí)，僅對(duì)包含UMI的轉(zhuǎn)錄本末端 (通常為3'末端)進(jìn)行性測序。

比對(duì)UMI條形碼

由于UMI數(shù)量（, N是UMIs的長度值）比每個(gè)細(xì)胞中的RNA分子數(shù)（~）少得多，每個(gè)UMI條形碼可能會(huì)連接到多個(gè)轉(zhuǎn)錄本，因此需要借助條形碼序列和reads比對(duì)位置兩個(gè)條件鑒定起始的轉(zhuǎn)錄本分子。第一步是比對(duì)UMI reads，推薦用STAR來處理，因?yàn)樗幚硭俣瓤烨逸敵龈哔|(zhì)量的BAM比對(duì)。此外，比對(duì)位置的準(zhǔn)確性對(duì)識(shí)別新的3'UTR區(qū)域很有意義。

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

處理UMI實(shí)驗(yàn)中的reads，通常有以下慣例：

UMI被添加到另一個(gè)配對(duì)read的序列名稱中。

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

Counting 條形碼

理論上，每個(gè)唯一的UMI-轉(zhuǎn)錄本對(duì)應(yīng)該對(duì)應(yīng)來源于一個(gè)RNA分子的所有reads。但是現(xiàn)實(shí)往往并非如此，最常見的原因是：

不同的UMI序列不一定表示它們是不同的UMI分子由于PCR或測序錯(cuò)誤，堿基替換可能導(dǎo)致新的UMI序列。較長的UMI出現(xiàn)堿基替換的機(jī)會(huì)更多。根據(jù)細(xì)胞條碼測序錯(cuò)誤估計(jì)，7-10％的10 bp長度的UMI中至少有一個(gè)堿基替換。如果錯(cuò)誤沒有糾正，將會(huì)過高估計(jì)轉(zhuǎn)錄本的數(shù)目。

不同的轉(zhuǎn)錄本不一定是不同的分子比對(duì)錯(cuò)誤或多個(gè)比對(duì)位置可能導(dǎo)致某些UMI對(duì)應(yīng)到錯(cuò)誤的基因/轉(zhuǎn)錄本。這種類型的錯(cuò)誤也會(huì)導(dǎo)致過高估計(jì)轉(zhuǎn)錄本的數(shù)目。

相同的UMI不一定意味著相同的分子UMI頻次的不同和短UMI可導(dǎo)致同一UMI和相同基因的不同mRNA分子相連，進(jìn)而可能低估轉(zhuǎn)錄本數(shù)量。

錯(cuò)誤糾正

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

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

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

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

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

這個(gè)方法的一個(gè)重要缺陷是它假設(shè)所有UMI出現(xiàn)頻率相同。但因?yàn)樾蛄蠫C含量不同引入的偏差使得這一假設(shè)在大多數(shù)情況下這是不正確的。

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

• UMI-tools

• PoissonUMIs

• zUMIs

• dropEst

下游分析

當(dāng)前的UMI平臺(tái)（DropSeq，InDrop，ICell8）展現(xiàn)出從低到高變化很大的捕獲效率，如下圖所示。

這一高可變性可能會(huì)引入很強(qiáng)的偏差，需要在下游分析時(shí)考慮到。現(xiàn)在的分析通常根據(jù)細(xì)胞類型或生物通路把細(xì)胞/gene混合一起增加檢測能力。更合適的統(tǒng)計(jì)分析方法亟待研究以便更好地調(diào)整這些偏差，使得結(jié)果更能反映真實(shí)現(xiàn)象。

練習(xí)1?數(shù)據(jù)是三個(gè)不同來源的誘導(dǎo)多功能干細(xì)胞的UMI counts和read counts?(有關(guān)此數(shù)據(jù)集的詳細(xì)信息請(qǐng)參閱后續(xù)文章)。

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

使用此數(shù)據(jù)：

繪制捕獲效率的變化

確定擴(kuò)增率：每個(gè)UMI對(duì)應(yīng)的平均reads數(shù)。

# 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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。