本文來引自倪帥科學(xué)網(wǎng)博客 http://blog.sciencenet.cn/blog-635619-884213.html 

人的基因組一共有兩萬多個(gè)基因，但是這些基因不是每時(shí)每刻都在表達(dá)，在不同發(fā)育時(shí)期和不同組織中，基因的表達(dá)是不同的，一個(gè)檢測這些表達(dá)的有效的方法就是RNA-seq，它結(jié)合了下一代測序的技術(shù)來對(duì)細(xì)胞整個(gè)的mRNA進(jìn)行測序，從而確定每一個(gè)基因的表達(dá)量和表達(dá)區(qū)段，主要用在分析不同條件下細(xì)胞內(nèi)基因表達(dá)差異和分析基因表達(dá)的不同可變剪接上。

RNAseq分析大致分下面幾個(gè)步驟，首先要把測到的序列map到基因組上，然后根據(jù)map到的區(qū)段對(duì)細(xì)胞構(gòu)建轉(zhuǎn)錄本，然后比較幾種細(xì)胞的轉(zhuǎn)錄本并且合并，最后衡量差異和可變間接和其他的分析。

1. Mapping

所有的序列分析的第一步，都是把測到的序列map到基因組上，這樣就能知道序列原來是在基因組的什么地方。mapping一般基于兩種快速索引算法，一種是哈希，MOSAIK，SOAP，SHRiMP用的就是這種算法，在對(duì)參照基因組建好哈希表之后，可以在常數(shù)次的運(yùn)算里查找到給定序列的位置，雖然高效，但是由于基因組有些區(qū)段重復(fù)性很高，所以查找次數(shù)雖是常數(shù)，但有時(shí)會(huì)變得非常大，降低效率；還有一種叫Burrows-Wheeler變換，BWA，Bowtie 和SOAP2都是用它，Burrows-Wheeler變換的設(shè)計(jì)比哈希更加巧妙，它最開始是一種文本壓縮算法，文本重復(fù)性越高，它的壓縮比就越大，這正好克服了基因組重復(fù)性高的問題，而且對(duì)于一個(gè)精確的序列查找，最多在給定序列的長度的次數(shù)里就能找到匹配，所以說基于Burrows-Wheeler變換的軟件在mapping里用得更加廣泛。可是RNAseq的map還有另外一個(gè)問題，那就是要允許可變剪接的存在，因?yàn)橐粭lRNA不一定是一個(gè)外顯子表達(dá)出來的，也有可能是幾個(gè)外顯子結(jié)合在了一起，原來基因里的內(nèi)含子被空了出來，這些內(nèi)含子的長度從五十到十萬個(gè)堿基不等，如果直接用DNAseq的方法的話去在基因組里尋找，有些正好在兩個(gè)exon連接處的序列就會(huì)有錯(cuò)配，而且有些在進(jìn)化過程中遺漏下來的假基因是沒有intron的，這樣就導(dǎo)致有些序列會(huì)被map到假基因上去，使假基因的表達(dá)變得很高，所以，傳統(tǒng)的bwa和bowtie在RNAseq里都不是最好的選擇。

更加適合RNA mapping的軟件需要克服上面的兩個(gè)問題，Tophat，subread, STAR, GSNAP, RUM, MapSplice都是為RNA測序而開發(fā)的，我只用過Tophat，它的新版Tophat2，在map的過程中分三個(gè)步驟，如果基因注釋文件存在的話，它會(huì)先用注釋文件的轉(zhuǎn)錄組來map，然后再對(duì)剩下的序列用bowtie進(jìn)行普通的map，最后再用bowtie里用過的所有的序列做剪接map，所以跟其他的軟件比起來會(huì)有比較高的正確率。在運(yùn)行tophat之前，要對(duì)參考基因組作index，samtools可以輕松搞定。

samtools faidx hg38.fasta

如果用的是tophat的話，還需要用bowtie2做index.

bowtie2-build genome.fa genome

hg38.fasta是人類的參考基因組，對(duì)參考基因組做index是為了提高mapping時(shí)查找的效率.

(對(duì)懶人來說，在這里可以直接下載到Bowtie打包做好的index文件，這樣就可以省略掉做index的步驟了 :D )

然后用Tophat開始mapping：

path/tophat

-p 8

-G $path_ref/Homo.GTF  

-o tophat_output

$path_ref/hg38

single_end.fastq;

-p指定用幾個(gè)線程來工作，-G指定注釋文件的位置，-o指定輸出文件的路徑和文件名，最后兩個(gè)參數(shù)分別告訴參考基因組的位置和要map的fastq文件。在參考基因組里不用加.fa的后綴，因?yàn)槌绦蜻€要去尋找其他后綴的index文件。

這步的運(yùn)行時(shí)間是根據(jù)fastq文件的大小和設(shè)定的線程數(shù)來決定的，一般單端8個(gè)線程需要每G一個(gè)小時(shí)，雙端各4G，線程數(shù)設(shè)為16的話需要五六小時(shí)，運(yùn)行完之后會(huì)在fastq文件的目錄里產(chǎn)生一個(gè)-o命令指定的文件夾，這個(gè)文件夾里有幾個(gè)bam文件和bed文件，還有一個(gè)summary，在下一步需要用到的是accepted_hits.bam這個(gè)文件。

注釋文件的后綴是GTF，它包含了所有已知的基因的外顯子在基因組中的位置。（hg38的注釋文件可以在這里下載）所以對(duì)于已經(jīng)map在基因組上的序列，我們可以直接根據(jù)它的位置從注釋文件里查找它是不是屬于一個(gè)外顯子，或者是一個(gè)轉(zhuǎn)錄本。對(duì)于要不要在這一步提供注釋文件，各有各的看法，我用單端測序的序列用兩種方法做了實(shí)驗(yàn)，發(fā)現(xiàn)他們有些差別：

有注釋的時(shí)候：

沒注釋的時(shí)候：

發(fā)現(xiàn)沒注釋的時(shí)候有更多的序列被map上去，但是重復(fù)的map也變多了，但是既然Tophat2三步map的第一步是根據(jù)注釋文件來map，我還是覺得在運(yùn)行tophat的時(shí)候用注釋文件比較好。

Mapping的最后一步是去除map到基因組中多于一處的序列，如果出現(xiàn)好幾個(gè)序列都map在完全相同的一個(gè)區(qū)段，那么就應(yīng)該只保留一個(gè)這樣的序列，所以，只保留匹配最高的那一個(gè)。而且這樣的序列占很大一部分，這步也很簡單，samtools里的rmdup可以輕松解決：

samtools rmdup -s input.bam output.bam

-s小寫是告訴samtools，bam文件是單端測序的結(jié)果，不指定-s的話默認(rèn)是雙端。

2。構(gòu)建轉(zhuǎn)錄本

Mapping完了以后，cufflinks就可以把map到基因組里的序列組裝成一個(gè)轉(zhuǎn)錄組了，這個(gè)轉(zhuǎn)錄組理論上包含了所有當(dāng)時(shí)細(xì)胞里的所有mRNA，組裝好的轉(zhuǎn)錄組包含了可能的剪切信息和所有轉(zhuǎn)錄的表達(dá)量，這個(gè)表達(dá)量是根據(jù)map到基因組的序列的總數(shù)和每個(gè)轉(zhuǎn)錄片斷的長度進(jìn)行歸一化的，聽起來比較難懂，它是對(duì)于在轉(zhuǎn)錄片斷里的每一千個(gè)堿基對(duì)，在每一百萬個(gè)成功map的序列中，map在這一千個(gè)堿基對(duì)上的序列的比例，fragments per kilobase of transcript per million mapped fragments (FKPM)。

FPKM是這么算出來的：

在公式里，C代表的是map在這一千個(gè)堿基對(duì)上的序列的個(gè)數(shù)，N是所有成功map的序列的個(gè)數(shù)，L是轉(zhuǎn)錄片斷的長度。

命令：

~/software/cufflinks-2.2.1/cufflinks

-g  /wrk/ref/Homo_sapiens.GRCh38.78.gtf

-o ../cufflinks_sample1

-p 8

accepted_hits.bam;

然后在輸出的cufflinks_sample1文件夾里會(huì)產(chǎn)生四個(gè)文件，genes.fpkm_tracking， isoforms.fpkm_tracking，skipped.gtf 和 transcripts.gtf，下一步需要用到的就是transcripts.gtf這個(gè)文件，transcripts.gtf就是這個(gè)樣品的轉(zhuǎn)錄組。

3。合并轉(zhuǎn)錄組

為了比較不同樣本間的差異，需要把實(shí)驗(yàn)組和對(duì)照組的轉(zhuǎn)錄組合并起來，cuffmerge不僅可以用來合并兩個(gè)或者多個(gè)轉(zhuǎn)錄組，還能把注釋過后的基因組的信息也合并起來，從而找到新的基因可變剪接，提高合并轉(zhuǎn)錄組的質(zhì)量。有人說需要在合并之前用cuffcompare，但從官網(wǎng)的說法來看沒有是必要的。它們最大的區(qū)別是，cuffcompare不改變?cè)袠颖纠锏霓D(zhuǎn)錄片段，只是將他們的位置作比較，輸出的combined文件也只是包含了所有的小轉(zhuǎn)錄片段，而cuffmerge會(huì)尋找?guī)讉€(gè)樣本間的不同，試著把幾個(gè)樣本里的轉(zhuǎn)錄片段從頭開始盡可能拼接成更長更完整的片段，所以cuffmerge的輸出merge文件比cuffcompare輸出的combined文件更有說服力。兩個(gè)小工具的作用都是為了生成一個(gè)合并的注釋文件給接下來要用的cuffdiff。

命令：

~/software/cufflinks-2.2.1/cuffmerge

-g /wrk/ref/Homo_sapiens.GRCh38.78.gtf

-s /wrk/ref/hg38.fa

-p 16

assembly_list.txt

在最后輸入的assemby_list.txt文件里，要寫上所有需要合并的transcripts.gtf文件的路徑，兩個(gè)和多個(gè)都可以，然后cuffmerge就會(huì)生成一個(gè)merged文件夾，里面有一個(gè)merged.gtf，這個(gè)文件就是合并好的轉(zhuǎn)錄組。

4. 基因差異表達(dá)分析

最后一步就是分析可變剪接和差異表達(dá)了，用到的小工具叫cuffdiff，這個(gè)的輸入比較復(fù)雜，不僅需要上一步的merge文件，還需要每個(gè)樣本的mapping結(jié)果的bam文件，最后還需要對(duì)每一個(gè)bam文件對(duì)應(yīng)的樣本按順序起一個(gè)名字作為標(biāo)簽，標(biāo)簽之間記得用逗號(hào)。

~/software/cufflinks-2.2.1/cuffdiff

-o diff22rv1/

-labels vap,vcap_neg_control

-p 32

-u merged_asm/merged.gtf  

vcap/accepted_hits.bam vcapneg/accepted_hits.bam

在這個(gè)例子里只有一個(gè)case和一個(gè)control，所以我們只要兩個(gè)標(biāo)簽，在最后按順序輸入bam文件。

diff的輸出比較多，他會(huì)對(duì)每個(gè)基因，每個(gè)轉(zhuǎn)錄片段，每個(gè)編碼序列，和每個(gè)基因的不同剪接體進(jìn)行FPKM，個(gè)數(shù)和樣本間差異進(jìn)行分析，最后生成幾組不同的文件，按照不同的分析需求，就可以試著往下分析了。

到現(xiàn)在結(jié)果就基本都差不多了，剩下的主要是作圖了，發(fā)現(xiàn)新的基因的可變剪接，發(fā)現(xiàn)差異表達(dá)的基因，對(duì)差異表達(dá)的基因做富集分析等等。作圖也是非常重要的環(huán)節(jié)，再好的結(jié)果也需要有好的圖表示出來，可變剪接我也沒做過，但是如果做差異表達(dá)的話，CummRbund是一個(gè)非常兼容cufflinks的作圖工具。

CummRbund是R里的一個(gè)包，用來分析cuffdiff的結(jié)果非常方便，在安裝好這個(gè)包之后，要做的只是把路徑改在cuffdiff生成結(jié)果的文件夾里，然后在R里運(yùn)行這兩行代碼就好了。

library(cummeRbud)      #這一行需要在R里加載已經(jīng)安裝好的cummeRbud包

cuff <- readCufflinks()    #這一行就告訴R把所有cuffdiff生成的結(jié)果讀入cuff這個(gè)變量里。

下一步就可以作差異表達(dá)的基因的熱圖了，這里稍微復(fù)雜一點(diǎn)：

#取前100個(gè)差異最顯著的基因，或者取多少個(gè)隨你便，標(biāo)準(zhǔn)是t檢驗(yàn)的p值，p值越小差異就越大。

gene.diff <- diffData(genes(cuff))

gene.diff.top <- gene.diff[order(gene.diff$p_value),][1:100,]

#找到這前100個(gè)差異基因的ID

myGeneIds <- gene.diff.top$gene_id

# 然后根據(jù)基因ID來得到基因的名字

myGenes <- getGenes(cuff, myGeneIds)

# 是畫圖

csHeatmap(myGenes, cluster="both")

作出的圖基本上是這個(gè)樣子的：

接下來就可以進(jìn)行富集分析了，有很多種方法，可以直接把基因的名字導(dǎo)出來后上傳到david來分析，也可以用bioconductor里的Goseq包來分析，詳細(xì)的goseq用法的代碼有點(diǎn)長，我以后會(huì)寫，如果需要的話可以直接去這里下載我使用Goseq時(shí)的代碼。

Cuff系列的分析流程是這篇文章介紹的，它里面也有非常詳細(xì)的命令和例子，可以去這里看看。Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks