RNA-seq數(shù)據(jù)分析

1、RNA-Seq數(shù)據(jù)分析 從原始的數(shù)據(jù)開始，進(jìn)行reads回帖，到拼接轉(zhuǎn)錄本，計算表達(dá)量，分析差異表達(dá)，最后可視化分析結(jié)果。     TopHat是一 個把reads回帖到基因組上的工具。首先用Bowtie把reads回帖到基因組上，然后通過拼接，我們就可以在基因組上看到一些reads堆疊起來的 區(qū)域，稱為consensus,這些consensus可能是一個真的外顯子，也有可能是幾個外顯子拼在一起的，或者一些別的情況。我們知道，經(jīng)典的剪切 位點一般都有GT和AG這樣的序列標(biāo)志，在consensus的邊界和內(nèi)部，TopHat會去找這樣的剪切位點，并且得到他們可能的組合。然后對于

2、那些沒 有被Bowtie貼到基因組上的reads,TopHat會對他們建立索引，去和這些可能的剪切位點比對，這樣就把跨越剪切位點的reads準(zhǔn)確地貼到基 因組上。    一些比較重要的命令行選項。       關(guān)于插入片段長度的選項：在RNA-Seq中，會把mRNA打斷成小的片段，然后對片段長度進(jìn)行iding篩選后拿去測序，如果選擇的片段長度是300bp，兩端各測序75bp的reads,中間的插入片段長度就應(yīng)該設(shè)為150bp.    下面是設(shè)置插入片段長度的標(biāo)準(zhǔn)差，如果選擇的片段長度比較集中，這個值可以設(shè)置的小一些

3、，反之應(yīng)該設(shè)置得大一些。    -G選項是提供哦呢一個已有的注釋文件。如果你分析的基因組被注釋得比較好了，最好能夠提供這個文件，這時TopHat就會先把reads往轉(zhuǎn)錄組上貼，沒有貼到轉(zhuǎn)錄組上的再往基因組上貼，最后把結(jié)果合并起來。我們知道大多數(shù)的轉(zhuǎn)錄組都是比基因組小得多的，而且junction reads可以直接貼到轉(zhuǎn)錄本上，所以這樣回帖的效力和準(zhǔn)確度都可以得到提高。       標(biāo)準(zhǔn)的Illumina平臺是不分鏈的，我們無法知道配對的reads哪個方向和轉(zhuǎn)錄本一致，哪個和轉(zhuǎn)錄本反向互補。對于分鏈的數(shù)據(jù)，也有兩 種情況，在f

4、irststrand這種分鏈方法中，第二個read和轉(zhuǎn)錄本方向一致，第一個read和轉(zhuǎn)錄本反向互補，在另一種fr- secondstrand分鏈方法中，就剛好反過來了。所以在分析的時候一定要弄清楚自己的數(shù)據(jù)有沒有分鏈，是怎么分鏈的。    下面是一個模擬的RNA-Seq數(shù)據(jù)集，雙端測序，有兩種處理，每種處理有3個重復(fù)，這里C代表處理，R代表重復(fù)，下面用C1R1進(jìn)行演示    首先，要有參考序列fasta文件，也就通常說的基因組序列。 TopHat是利用Bowtie2回帖reads,我們首先需要建立Bowtie2的索引文件： bowtie2-bu

5、ild genome.fa(基因組文件) genome  （注意程序和文件所在目錄）    我們還需要reads的fastq文件，雙端測序的數(shù)據(jù)，兩個fastq文件分別以下劃線1和2這樣的形式結(jié)尾。在實際分析中，需要對拿到的數(shù)據(jù)進(jìn)行質(zhì)量 評估和過濾等依稀類預(yù)處理工作，這些工作都是非常重要的。需要準(zhǔn)備注釋文件，當(dāng)然它不是必須的。它可以是GTF或者GFF3格式的文件，對于注釋得比較好 的基因組，在UCSC可以下載。   準(zhǔn)備好后就可以運行tophat了，-p是線程數(shù)，-G是注釋文件，-o是輸出文件夾，選項之后就是參考序列的索引，最后是兩個reads

6、的fastq文件。   看里面生成的文件，align-summary文件愛你，這個文件是reads回帖的一些統(tǒng)計信息。90%以上的回帖比例就非常好了，當(dāng)然百 分之70以上一個可以接受的范圍。bam文件詳細(xì)記錄了reads回帖到基因組上的情況，由于這是一個二進(jìn)制的文件，我們需要用samtools查看它。    Cufflinks是一套拼接轉(zhuǎn)錄本，計算表達(dá)量，計算差異表達(dá)的工具。盡可能拼接處最優(yōu)可能的轉(zhuǎn)錄本的結(jié)構(gòu)，并且估計它的表達(dá)量。 -G是提供一個注釋文件，并且告訴Cufflinks不要去拼接新的轉(zhuǎn)錄本，只能用注釋文件里提供的轉(zhuǎn)錄本。 

7、-g 也是提供一個注釋文件，但是Cufflinks會在這些已知轉(zhuǎn)錄本的指導(dǎo)下，拼接新的轉(zhuǎn)錄本。-u是告訴Cufflinks用更準(zhǔn)確的方法去處理貼到多個位點上的reads，如果沒有-u,Cufflinks只會把這些reads簡單地平均分配。 比如一個read貼到了10個位置，那么每個位置分得十分之一。加-u后會先進(jìn)行平均分配，然后按照這10個位置各自的表達(dá)量，計算read被分配到每個 位置的概率。實際上Cufflinks會用EM算法進(jìn)行迭代，計算在觀察到當(dāng)前數(shù)據(jù)的情況下，最優(yōu)可能的reads分配。library type和TopHat里面差不多。這里的bam文件就是剛才TopHat運行的結(jié)果。C

8、uffmerge當(dāng)我們使用Cufflinks處理多個數(shù)據(jù)之后，我們需要將其轉(zhuǎn)錄本數(shù)據(jù)整合為一個全面的轉(zhuǎn)錄本集合，Cuffmerge是一個將Cufflinks生成 的gtf文件融合為一個更加全面的轉(zhuǎn)錄本注釋結(jié)果的工具。如下圖所示，圖中的6個轉(zhuǎn)錄本被整合為一個轉(zhuǎn)錄本集合。同時我們可以利用基因組注釋文件，獲得更 加準(zhǔn)確可靠的結(jié)果。合并后的轉(zhuǎn)錄本集合為計算每個基因和轉(zhuǎn)錄本的表達(dá)量提供了一個統(tǒng)一的基礎(chǔ)。-g參數(shù) 指向參考GTF文件。   -p參數(shù)決定線程數(shù)   -s參考指向基因組DNA序列。如果是一個文件夾，每個contig是一個fasta文件。如果是一個fast

9、a文件，則所有的contigs都需要在里面。    最后一項是一個列表，內(nèi)容包括經(jīng)過Cufflinks拼接的轉(zhuǎn)錄本的文件路徑。流程如下：首先，我們需要用cat命令創(chuàng)建一個所有拼接的轉(zhuǎn)錄本的文件路徑列表。  cat> assemblies.txt  復(fù)制上圖的列表  然后運行cuffmerge,運行后的結(jié)果存儲在merge_asm這個文件夾里面。其文件夾內(nèi)包含一個Logs文件夾以及一個.gtf文件，也就是我們經(jīng)過整合的轉(zhuǎn)錄本文件。cuffdiff    當(dāng)我們利用Cufflinks獲得了拼接的轉(zhuǎn)錄本時，我們就可以計算不同

10、樣品中轉(zhuǎn)錄本的表達(dá)量。計算的簡單原理在于測序深度和外顯子長度一定時，Read 的數(shù)量與對應(yīng)的轉(zhuǎn)錄本數(shù)量成正比。通過對Reads進(jìn)行計數(shù)計算轉(zhuǎn)錄本的表達(dá)量。同時cuffdiff可以計算不同條件下轉(zhuǎn)錄本表達(dá)水平的顯著性差異。    -u命令指cuffdiff對回帖的基因組中多個位置的read進(jìn)行一個初步的估計，然后加權(quán)分配到各個基因組位置。而不是簡單的平均分配，其功能與Cufflinks中的u命令相同。   -L為每個樣品標(biāo)上名稱。    接下來Cuffmerge產(chǎn)生的gtf文件，Cuffdiff需要它提供的注釋進(jìn)行初始轉(zhuǎn)錄產(chǎn)物和可變剪切等定量分析。最后是TopHat產(chǎn)生的bam文件，如果一個樣品中有多個實驗重復(fù)，那么我們需要提供哦呢bam文件列表，文件名之間以逗號隔開。        運行之后，cuffdiff輸出的文件在diff_out目錄之下。其中包括一些按類別統(tǒng)計的表達(dá)水平結(jié)果，如果有相同的轉(zhuǎn)錄起始位點，或具有相同的編碼區(qū)的轉(zhuǎn)錄本的表達(dá)水平，我么你可以利用他們進(jìn)行下一步的分析。CummeRbund當(dāng)我們對