原文來(lái)自：聊聊轉(zhuǎn)錄組測(cè)序——2.數(shù)據(jù)分析與解讀（上）
流程概覽
轉(zhuǎn)錄組測(cè)序的分析流程大致可以分成三類逢渔，包括基因組比對(duì)（Genome mapping）胃夏、轉(zhuǎn)錄組比對(duì)（Transcriptome mapping）、轉(zhuǎn)錄組組裝（Reference-free assembly）续镇，見下圖。其中第三種主要是用于分析沒(méi)有參考基因組和基因注釋的物種，應(yīng)用場(chǎng)合較少且不適合新手入門梦碗。對(duì)于人、小鼠蓖救、大鼠等模式物種洪规，通常用前兩種方法進(jìn)行分析。雖然轉(zhuǎn)錄組比對(duì)相關(guān)軟件和流程同樣層出不窮循捺，但對(duì)于基因組信息較為完善的模式物種斩例，推薦使用基因組比對(duì)的方式進(jìn)行分析，具體原因下文的“比對(duì)”部分會(huì)有說(shuō)明从橘。我們下面也主要對(duì)基因組比對(duì)的方法進(jìn)行介紹念赶。
圖片來(lái)源：https://www.ncbi.nlm.nih.gov/pubmed/26813401
1. Data Cleaning
從原始數(shù)據(jù)（Raw Data）到干凈數(shù)據(jù)（Clean
Data）的過(guò)程，有人翻譯成“數(shù)據(jù)清洗”恰力，實(shí)在叫不習(xí)慣叉谜，那我就不翻譯了。
Illumina測(cè)序儀下機(jī)的數(shù)據(jù)通常為Bcl格式踩萎，是將同一個(gè)測(cè)序通道（Lane）所有樣品的數(shù)據(jù)混雜在一起的停局，所以公司一般不會(huì)提供Bcl文件。測(cè)序公司使用Illumina官方出品的Bcl2FastQ軟件香府，根據(jù)Index序列分割轉(zhuǎn)換成每個(gè)樣品的FastQ文件董栽，打開長(zhǎng)這樣：
每一條序列（read）包含四行，第一行是read的ID企孩，第二行是序列锭碳，第四行是序列中每個(gè)堿基的測(cè)序質(zhì)量（更具體的細(xì)節(jié)可參考FASTQ format - Wikipedia）。
原始數(shù)據(jù)沒(méi)法直接分析柠硕，是因?yàn)椴糠謗eads測(cè)序質(zhì)量較低工禾，可能會(huì)誤導(dǎo)后續(xù)結(jié)果运提，因此需要對(duì)低質(zhì)量堿基太多或N（未能識(shí)別的堿基）太多的reads進(jìn)行去除；此外闻葵，部分測(cè)序文庫(kù)的插入片段太短民泵，導(dǎo)致測(cè)到兩側(cè)的接頭序列（請(qǐng)參考上一篇的測(cè)序文庫(kù)結(jié)構(gòu)圖理解），這些序列接頭也需要從reads中去除槽畔。最后栈妆，我們也會(huì)對(duì)清洗前后的Raw Data和Clean Data進(jìn)行評(píng)估，評(píng)估內(nèi)容包括堿基質(zhì)量厢钧、序列長(zhǎng)度鳞尔、堿基比例、GC含量早直、重復(fù)序列寥假、Kmers等（詳情請(qǐng)參考FastQC說(shuō)明文檔FastQC A Quality Control tool for High Throughput Sequence Data）。
最后說(shuō)一句霞扬，其實(shí)大多數(shù)測(cè)序公司都會(huì)提供CleanData的糕韧。
#常用軟件#
我以前都是用cutadapt + FASTX-Toolkit的組合，直到同事們給我推薦了Trim Galore喻圃，質(zhì)量評(píng)估使用FastQC萤彩。
2. 比對(duì)
由于二代測(cè)序的reads長(zhǎng)度通常介于50~300個(gè)堿基，因此即便使用雙端測(cè)序斧拍，也基本不可能覆蓋完整的mRNA轉(zhuǎn)錄本雀扶，因此想直接用FastQ文件從頭分析測(cè)到了哪些轉(zhuǎn)錄本需要非常復(fù)雜的分析和計(jì)算。好在通常情況下肆汹，公共數(shù)據(jù)庫(kù)已經(jīng)提供了測(cè)序樣品的基因組和轉(zhuǎn)錄本的序列愚墓。因此我們只需要知道，每一條reads來(lái)自哪一條轉(zhuǎn)錄本就可以了昂勉，這個(gè)將reads與參考（Reference）基因組/轉(zhuǎn)錄組的序列進(jìn)行比較和匹配的過(guò)程转绷，我們通常稱之為“比對(duì)”（文獻(xiàn)中提到的read alignment和mapping通常說(shuō)的都是這個(gè)）。
正如前文所述硼啤，轉(zhuǎn)錄組測(cè)序的比對(duì)通常分為基因組比對(duì)和轉(zhuǎn)錄組比對(duì)兩種，顧名思義斧账，基因組比對(duì)就是把reads比對(duì)到完整的基因組序列上谴返，而轉(zhuǎn)錄組比對(duì)則是把reads比對(duì)到所有已知的轉(zhuǎn)錄本序列上。如果不是很急或者只想知道已知轉(zhuǎn)錄本表達(dá)量咧织，個(gè)人建議使用基因組比對(duì)的方法進(jìn)行分析嗓袱，理由如下：
① 轉(zhuǎn)錄組比對(duì)需要準(zhǔn)確的已知轉(zhuǎn)錄本的序列，對(duì)于來(lái)自未知轉(zhuǎn)錄本（比如一些未被數(shù)據(jù)庫(kù)收錄的lncRNA）或序列不準(zhǔn)確的reads無(wú)法正確比對(duì)习绢；
② 與上一條類似渠抹，轉(zhuǎn)錄組比對(duì)不能對(duì)轉(zhuǎn)錄本的可變剪接進(jìn)行分析蝙昙，數(shù)據(jù)庫(kù)中未收錄的剪接位點(diǎn)會(huì)被直接丟棄；
③ 由于同一個(gè)基因存在不同的轉(zhuǎn)錄本梧却，因此很多reads可以同時(shí)完美比對(duì)到多個(gè)轉(zhuǎn)錄本奇颠，reads的比對(duì)評(píng)分會(huì)偏低，可能被后續(xù)計(jì)算表達(dá)量的軟件舍棄放航，影響后續(xù)分析（有部分軟件解決了這個(gè)問(wèn)題）烈拒；
④ 由于與DNA測(cè)序使用的參考序列不同，因此不利于RNA和DNA數(shù)據(jù)的整合分析广鳍。
而上面的問(wèn)題使用基因組比對(duì)都可以解決荆几。
此外，值得注意的是赊时，RNA測(cè)序并不能直接使用DNA測(cè)序常用的BWA吨铸、Bowtie等比對(duì)軟件，這是由于真核生物內(nèi)含子的存在祖秒，導(dǎo)致測(cè)到的reads并不與基因組序列完全一致（如下圖所示）诞吱，因此需要使用Tophat/HISAT/STAR等專門為RNA測(cè)序設(shè)計(jì)的軟件進(jìn)行比對(duì)。
圖片來(lái)源：GATK | Best Practices
比對(duì)結(jié)果會(huì)展示為BAM/SAM文件狈涮，其中BAM格式是SAM格式的二進(jìn)制版本（請(qǐng)理解為壓縮后的版本狐胎，用Samtools可以打開），打開之后長(zhǎng)這樣：
BAM文件中每行代表一條reads的比對(duì)信息歌馍，其中第一列是read的ID握巢，第二列為FLAG（包括是否雙端比對(duì)，比對(duì)位點(diǎn)是否唯一等信息）松却，第三列為比對(duì)的染色體暴浦，第四列為比對(duì)的起始位置，第六列為CIGAR值晓锻，代表比對(duì)的具體方式（例60M2D80M代表60個(gè)堿基完美匹配+2個(gè)堿基缺失+80個(gè)堿基完美匹配）等等歌焦，BAM文件的具體內(nèi)容可參考SAM - Genome Analysis Wiki和http://samtools.github.io/hts-specs/SAMv1.pdf（后面這個(gè)要翻墻）。
#常用軟件#
基因組比對(duì)：
Tophat2：可以說(shuō)是最被公認(rèn)的RNA測(cè)序比對(duì)軟件（實(shí)際上是在DNA比對(duì)軟件Bowtie的基礎(chǔ)上做了一個(gè)殼）砚哆，相信很多做RNA測(cè)序的同學(xué)都是看著Tophat發(fā)表在Nature Protocol上的步驟一步步入門RNA測(cè)序的独撇；
HISAT2：Tophat2的非正式升級(jí)版本（因?yàn)閾?jù)說(shuō)還會(huì)有Tophat3），在Tophat的算法基礎(chǔ)了上做了大量的改進(jìn)躁锁，而且克服了Tophat最大的缺點(diǎn)——速度慢纷铣，Nature Protocol上同樣發(fā)表了操作流程；
STAR：ENCODE計(jì)劃御用比對(duì)軟件战转，權(quán)威程度可以與Tophat平起平坐搜立，并且比對(duì)速度極快；
MapSplice：TCGA使用的比對(duì)軟件槐秧，我自己沒(méi)用過(guò)啄踊；
RSEM：RSEM更像一個(gè)軟件包而不是一個(gè)比對(duì)軟件忧设，能夠提供從比對(duì)到計(jì)算差異表達(dá)的所有步驟，由于不需要自己寫代碼串聯(lián)不同軟件生成的數(shù)據(jù)格式颠通，因此用起來(lái)比較省時(shí)省力址晕，值得注意的是，TCGA使用MapSplice比對(duì)后再用RSEM計(jì)算表達(dá)量蒜哀，并沒(méi)有直接只用RSEM原裝的Bowtie的比對(duì)結(jié)果斩箫。
轉(zhuǎn)錄組比對(duì)：這類型的軟件我用的不多，最近嘗試過(guò)Nature Methods上面發(fā)表的Salmon撵儿，能從Clean reads直接算到表達(dá)量乘客，優(yōu)點(diǎn)是，快淀歇，非骋缀耍快。然而這個(gè)軟件連BAM文件都沒(méi)生成浪默，雖然只是定量的話BAM文件的確沒(méi)什么用就是了…
1.5 #可選步驟# 核糖體RNA（rRNA）去除
嗯牡直，寫完2再寫1.5是我不對(duì)。
如果對(duì)上一篇還有印象的話纳决，我們?cè)岬脚鲆荩D(zhuǎn)錄組測(cè)序有一種 偏貴的 使用核糖體RNA去除技術(shù)構(gòu)建文庫(kù)的測(cè)序。但是經(jīng)常做實(shí)驗(yàn)的你一定知道阔加，這種去除是沒(méi)法做到100%去除rRNA的饵史，更糟糕的是，同一批測(cè)序的每個(gè)樣品胜榔，rRNA的去除效率也會(huì)有一定差別的胳喷！
由于rRNA都是非編碼RNA序列，因此如果我們后續(xù)分析需要使用轉(zhuǎn)錄本組裝的方法鑒別新的lncRNA（long non-coding RNA夭织，長(zhǎng)非編碼RNA）吭露，這些rRNA的序列特征很容易對(duì)lncRNA的鑒定造成干擾，因此我們必須對(duì)這些rRNA序列進(jìn)行去除尊惰。
當(dāng)然讲竿，如果不涉及組裝新lncRNA的話，rRNA的存在對(duì)分析結(jié)果的影響并不大弄屡。但如果樣品間rRNA殘留率差別較大戴卜，對(duì)定量的準(zhǔn)確性會(huì)有較大影響，因此有能力的話還是建議去除rRNA序列琢岩。就算不去除，用比對(duì)軟件算算rRNA序列占總數(shù)據(jù)量的比例也是好的师脂，一旦不小心發(fā)現(xiàn)12G的數(shù)據(jù)里面6G都來(lái)自于rRNA…（嗯我不是教你們跟公司撕X…）
#常用軟件#
核糖體去除實(shí)際上也是通過(guò)比對(duì)來(lái)進(jìn)行担孔，我在Rfam上下載rRNA的序列后江锨，直接使用Bowtie2進(jìn)行比對(duì)。
至于比對(duì)核糖體之后怎么拿到?jīng)]有rRNA的FastQ文件糕篇，我不太清楚別人是怎么做的啄育，我是用Python把沒(méi)比對(duì)上的Reads的ID提取出來(lái)存成一個(gè)表格，再用Seqtk提取FastQ文件拌消。
#########################未完待續(xù)挑豌，我們下回分解#########################
最后關(guān)于常用軟件：
首先推薦OMIC TOOLS網(wǎng)站（https://omictools.com/rna-seq-category），上面收集了大量高通量測(cè)序相關(guān)的軟件以及軟件之間的對(duì)比評(píng)測(cè)文獻(xiàn)墩崩。
其次氓英，我寫的常用軟件基本基于我自己和身邊同事的使用情況，以及在文獻(xiàn)中看到的情況鹦筹，并不涉及對(duì)軟件自身性能的評(píng)判铝阐。比如華大基因開發(fā)的比對(duì)用軟件SOAPsplice，沒(méi)有提到不是說(shuō)它不好铐拐，而是我的確沒(méi)用過(guò)徘键，而且用的人也比較少。
至于如何判斷哪個(gè)軟件更好遍蟋，可以參考軟件評(píng)測(cè)的文獻(xiàn)吹害，但沒(méi)必要以此為絕對(duì)標(biāo)準(zhǔn)凉夯，例如我見到有評(píng)測(cè)說(shuō)STAR完爆Tophat的运悲，也有說(shuō)Tophat完爆STAR的，而且這兩篇評(píng)測(cè)都發(fā)在Nature Method上经柴。我個(gè)人認(rèn)為能一直存活下來(lái)也被廣泛使用的軟件大多各有千秋挟憔，對(duì)于一個(gè)沒(méi)精力看源代碼（也看不懂java和C++）的使用者來(lái)說(shuō)钟些，想要評(píng)估一個(gè)軟件的好壞，直接參考高分雜志的使用情況來(lái)推算绊谭，可能是更靠譜的選擇政恍。
