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分享
構(gòu)建表達(dá)矩陣
scRNA-seq數(shù)據(jù)的許多分析以表達(dá)矩陣為起點(diǎn)。一般來(lái)講温自，表達(dá)矩陣的每一行代表一個(gè)基因玄货，每一列代表一個(gè)細(xì)胞（但是一些作者會(huì)做個(gè)轉(zhuǎn)置）。每個(gè)條目代表特定基因在給定細(xì)胞中的表達(dá)水平悼泌。而表達(dá)值的測(cè)量單位取決于建庫(kù)方案和所用的標(biāo)準(zhǔn)化方法松捉。
reads質(zhì)控
見(jiàn)前面章節(jié)FastQC部分。
另外馆里，使用Integrative Genomics Browser（IGV）或SeqMonk通常對(duì)數(shù)據(jù)可視化很有幫助隘世，具體見(jiàn)下可柿。
Reads比對(duì)
見(jiàn)前面章節(jié)STAR部分和Kallisto部分。
注釋的很好的模式生物（例如小鼠丙者，人）有著大量全長(zhǎng)轉(zhuǎn)錄本數(shù)據(jù)集复斥，偽比對(duì)方法（例如Kallisto，Salmon)可能優(yōu)于常規(guī)比對(duì)方法械媒。drop-seq方法獲得的數(shù)據(jù)集有數(shù)以千萬(wàn)條reads目锭，偽比對(duì)工具的運(yùn)行時(shí)間比傳統(tǒng)比對(duì)工具會(huì)少幾個(gè)數(shù)量級(jí)，更有時(shí)間優(yōu)勢(shì)纷捞。從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序列合并作為共同”參考基因組”缀辩。具體見(jiàn)之前文件格式部分臭埋。
注意，使用UMI時(shí)臀玄，應(yīng)從read序列中刪除其條形碼瓢阴。常見(jiàn)的是將條形碼加到read名稱(chēng)上。
一旦reads完成了到基因組的比對(duì)健无，我們需要檢查比對(duì)率和確保有足夠多的reads比對(duì)回了參考基因組荣恐。根據(jù)我們的經(jīng)驗(yàn)，小鼠或人類(lèi)細(xì)胞中read的比對(duì)率為60-70％累贤。但是這個(gè)結(jié)果可能會(huì)因建庫(kù)方案叠穆、read長(zhǎng)度和比對(duì)工具參數(shù)設(shè)置而有所不同。常規(guī)上臼膏，我們希望所有細(xì)胞都具有相似的read比對(duì)率硼被。如果有樣品比對(duì)率異常低或比對(duì)回去的reads異常低，則需要多加注意甚至從后續(xù)分析中移除渗磅。較低的read比對(duì)率通常表示存在污染嚷硫。生信寶典建議取出幾十條未必對(duì)回去的reads做個(gè)blast，看下能否比對(duì)到其它物種來(lái)確定是污染還是測(cè)序錯(cuò)誤還是比對(duì)參數(shù)設(shè)置的問(wèn)題始鱼。
一個(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ì)胞測(cè)序中長(zhǎng)基因的表達(dá)做了過(guò)度校正医清，因此我們建議使用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ì)保留下來(lái)，但后期細(xì)胞質(zhì)控時(shí)這兩個(gè)細(xì)胞會(huì)因?yàn)楹颂求wRNA reads比例過(guò)高而移除残腌。
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開(kāi)發(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)估取決于采用的建庫(kù)方案，例如許多scRNA-seq用poly-A selection捕獲轉(zhuǎn)錄本讥巡。這個(gè)方法可以排除核糖體RNA污染掀亩，但會(huì)導(dǎo)致3'區(qū)域更容易測(cè)到。下圖展示了測(cè)序reads分布的3'偏好性欢顷，和去除的三個(gè)異常細(xì)胞的結(jié)果 (應(yīng)該是最下面3條槽棍，推測(cè)是降解嚴(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)錄過(guò)程中添加到轉(zhuǎn)錄本上的短的（4-10 bp）隨機(jī)條形碼序列碴犬。它們使得測(cè)序reads可以對(duì)應(yīng)到單個(gè)轉(zhuǎn)錄本题暖，從而去除擴(kuò)增噪聲和偏好性按傅。
當(dāng)測(cè)序含UMI的文庫(kù)時(shí)，僅對(duì)包含UMI的轉(zhuǎn)錄本末端 (通常為3'末端)進(jìn)行性測(cè)序芙委。
比對(duì)UMI條形碼
由于UMI數(shù)量（, N是UMIs的長(zhǎng)度值）比每個(gè)細(xì)胞中的RNA分子數(shù)（~）少得多逞敷，每個(gè)UMI條形碼可能會(huì)連接到多個(gè)轉(zhuǎn)錄本，因此需要借助條形碼序列和reads比對(duì)位置兩個(gè)條件鑒定起始的轉(zhuǎn)錄本分子灌侣。第一步是比對(duì)UMI reads推捐，推薦用STAR來(lái)處理，因?yàn)樗幚硭俣瓤烨逸敵龈哔|(zhì)量的BAM比對(duì)侧啼。此外牛柒，比對(duì)位置的準(zhǔn)確性對(duì)識(shí)別新的3'UTR區(qū)域很有意義堪簿。
UMI測(cè)序通常由雙端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的序列名稱(chēng)中滴须。
reads按細(xì)胞條形碼分類(lèi)到單獨(dú)的文件中 (見(jiàn)前面的文章)舌狗。但對(duì)于細(xì)胞量極大的低深度測(cè)序數(shù)據(jù)集 (drop-seq)，可以將細(xì)胞條形碼添加到read名稱(chēng)中而不是拆分為單獨(dú)文件以減少文件數(shù)量扔水。
Counting 條形碼
理論上痛侍，每個(gè)唯一的UMI-轉(zhuǎn)錄本對(duì)應(yīng)該對(duì)應(yīng)來(lái)源于一個(gè)RNA分子的所有reads。但是現(xiàn)實(shí)往往并非如此魔市，最常見(jiàn)的原因是：
不同的UMI序列不一定表示它們是不同的UMI分子由于PCR或測(cè)序錯(cuò)誤主届，堿基替換可能導(dǎo)致新的
UMI序列。較長(zhǎng)的UMI出現(xiàn)堿基替換的機(jī)會(huì)更多待德。根據(jù)細(xì)胞條碼測(cè)序錯(cuò)誤估計(jì)君丁，7-10％的10 bp長(zhǎng)度的UMI中至少有一個(gè)堿基替換。如果錯(cuò)誤沒(méi)有糾正将宪，將會(huì)過(guò)高估計(jì)轉(zhuǎn)錄本的數(shù)目谈截。不同的轉(zhuǎn)錄本不一定是不同的分子比對(duì)錯(cuò)誤或多個(gè)比對(duì)位置可能導(dǎo)致某些
UMI對(duì)應(yīng)到錯(cuò)誤的基因/轉(zhuǎn)錄本。這種類(lèi)型的錯(cuò)誤也會(huì)導(dǎo)致過(guò)高估計(jì)轉(zhuǎn)錄本的數(shù)目涧偷。相同的UMI不一定意味著相同的分子UMI頻次的不同和短UMI可導(dǎo)致同一UMI和相同基因的不同mRNA分子相連簸喂，進(jìn)而可能低估轉(zhuǎn)錄本數(shù)量。
錯(cuò)誤糾正
如何最好的校正UMIs中的這些問(wèn)題仍然是一個(gè)活躍的研究領(lǐng)域燎潮。我們自己認(rèn)為的最好的解決上述問(wèn)題的方法是：
UMI-tools’喻鳄，設(shè)計(jì)了
directional-adjacency算法，同時(shí)考慮錯(cuò)配數(shù)目和相似UMI的相對(duì)頻率來(lái)識(shí)別PCR和測(cè)序錯(cuò)誤确封。(alevin, cellranger都是不錯(cuò)的選擇除呵，后面詳細(xì)介紹)問(wèn)題還無(wú)法完全解決，但通過(guò)刪除只有很少read支持的
UMIs-轉(zhuǎn)錄本對(duì)爪喘，或者移除所有多比對(duì)位置的reads颜曾，可能會(huì)減輕該問(wèn)題。Simple saturation (也稱(chēng)為”collision probability”)方法來(lái)估計(jì)分子的數(shù)量
其中N=唯一UMI條形碼的總數(shù)秉剑，n=觀(guān)察到的條形碼數(shù)泛豪。
這個(gè)方法的一個(gè)重要缺陷是它假設(shè)所有UMI出現(xiàn)頻率相同。但因?yàn)樾蛄蠫C含量不同引入的偏差使得這一假設(shè)在大多數(shù)情況下這是不正確的。
如何最好地處理和使用UMI在目前生物信息學(xué)界是一個(gè)活躍的研究領(lǐng)域诡曙。而我們了解到的幾種最近開(kāi)發(fā)的方法有：
UMI-tools
PoissonUMIs
zUMIs
dropEst
下游分析
當(dāng)前的UMI平臺(tái)（DropSeq臀叙，InDrop，ICell8）展現(xiàn)出從低到高變化很大的捕獲效率价卤，如下圖所示劝萤。
這一高可變性可能會(huì)引入很強(qiáng)的偏差，需要在下游分析時(shí)考慮到∩麒担現(xiàn)在的分析通常根據(jù)細(xì)胞類(lèi)型或生物通路把細(xì)胞/gene混合一起增加檢測(cè)能力床嫌。更合適的統(tǒng)計(jì)分析方法亟待研究以便更好地調(diào)整這些偏差，使得結(jié)果更能反映真實(shí)現(xiàn)象胸私。
練習(xí)1 數(shù)據(jù)是三個(gè)不同來(lái)源的誘導(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
