RNAseq定量分析方案

1、福建農(nóng)林大學基因組與生物技術研究中心整理 2013-12-31RNAseq定量分析方案RNAseq定量分析方案1一、實驗目的：2二、實驗大致流程2三、實驗前的準備活動23.1準備數(shù)據(jù)23.2確定涉及軟件是否安裝完畢，及其輸入輸出文件。3四、實驗過程44.1.利用bowtie2-bulid命令根據(jù)提供的參考基因組序列建立對應的索引文件。54.2.利用tophat命令將分別將代比較的reads maping到參考基因組序列上64.3利用cufflinks軟件分別將待測樣品的轉(zhuǎn)錄組reads拼接起來，并同時計算每個樣品各個基因的rpkm值74.4.利用coffmerge和cuffdiff軟件計算每個

2、樣品各個基因的fpkm值。8五利用R軟件查看結果文件。10一、實驗目的：利用已有的水稻基因組數(shù)據(jù)對來自兩棵不同的水稻進行轉(zhuǎn)錄組水平差異的研究。二、實驗大致流程1. 利用bowtie2-bulid命令根據(jù)提供的參考基因組序列建立對應的索引文件2. 利用tophat命令將分別將代比較的reads maping到參考基因組序列上3. 利用cufflinks軟件分別將待測樣品的轉(zhuǎn)錄組reads拼接起來，計算每個樣品各個基因的fpkm值4. 利用coffmerge和cuffdiff軟件計算每個樣品各個基因的fpkm值。5. 利用R軟件查看比較結果。三、實驗前的準備活動3.1準備數(shù)據(jù)像大多數(shù)生物實驗一樣，

3、做生物信息學實驗之前也是需要事先準備好“藥品”和“儀器”，不然到了關鍵時刻也是一樣會手忙腳亂的。現(xiàn)在我們先來談一談生物信息實驗需要準備的“藥品”數(shù)據(jù)。對于RNAseq實驗而言，這里至少需要準備一下幾個文件： 1.參考基因組序列文件 如 refrence.fa 參考基因組數(shù)據(jù)： 2.參考基因組注釋文件 如 refrence.gtf a參考基因組序列設為refrence.fa b參考基因注釋設為refrence.gtf F1_sample_R1.fastq待比較樣品F1：(為RNA測序數(shù)據(jù)) F1_sample_R1.fastqP1_sample_R1.fastq待比較樣品P2：(為RNA測序數(shù)據(jù)

4、) P1_sample_L007_R2.fastq3.2確定涉及軟件是否安裝完畢，及其輸入輸出文件。通過上面的大致流程，我們知道，我們至少需要用到以下幾個軟件，查看你的系統(tǒng)是否安裝了這些軟件只需要直接輸入這些命令就可以了：1.bowtie2-build http:/computing.bio.cam.ac.uk/local/doc/bowtie2.html Usage: bowtie2-build options* <reference_in> <bt2_index_base> reference_in comma-separated list of files wit

5、h ref sequences bt2_index_base write .bt2 data to files with this dir/basename2.tophat /tutorial.shtml#tophTopHat maps short sequences from spliced transcripts to whole genomes.Usage: tophat options <bowtie_index> <reads1,reads2,.> reads1,reads2,. quals1,quals2,.

6、 quals1,quals2,.3.cufflinks /manual.htmlUsage: cufflinks options <hits.sam>4.cuffmerge /manual.htmlUsage:cuffmerge Options <assembly_GTF_list.txt>5.cuffdiff /manual.htmlUsage: cuffdiff options <transcr

7、ipts.gtf> <sample1_hits.sam> <sample2_hits.sam> . sampleN_hits.samSupply replicate SAMs as comma separated lists for each condition: sample1_rep1.sam,sample1_rep2.sam,.sample1_repM.sam由于這些命令的參數(shù)都不少，所以不能一一在這里講解，想要了解這些軟件的參數(shù)作用最好的方法就是直接讀這些軟件的manual。四、實驗過程本實驗數(shù)據(jù)存儲在/share/Public/BioinfoTraim下

8、面的Train_RNAseq目錄里。如果你想做練習的話，只要在/share/Public/BioinfoTraim將數(shù)據(jù)重新拷貝一份，然后再拷貝的數(shù)據(jù)下進行操作就可以了。方法如下所示：現(xiàn)在正式進入RNAseq的實驗：1.進入到Train_RNAseq這個目錄里面(你則進入到自己拷貝的目錄里面)，查看一下該目錄下都有些什么文件。我們可以看到，在該目錄下，有我們需要的所有數(shù)據(jù)文件，說明數(shù)據(jù)已經(jīng)準備就緒。4.1.利用bowtie2-bulid命令根據(jù)提供的參考基因組序列建立對應的索引文件。我們知道參考基因組的序列文件時Osativa.fa文件，所以這就是bowite2-bulid命令的輸入文件，而b

9、owtie2-build的輸出文件時什么呢？根據(jù)bowtie2-build的Usage: bowtie2-build options* <reference_in> <bt2_index_base> reference_in comma-separated list of files with ref sequences bt2_index_base write .bt2 data to files with this dir/basename可知，bowtie2-build的輸出文件名可以由我們制定的，實際上是輸出文件不止一個，我們制定的只是這些輸出文件的文件名的前綴

10、而已。照此說來，我們這里可以用這句命令來實現(xiàn)我們的目標：# bowtie2-build Osativa.fa indexname但是你會看到輸出的軟件，把我的實驗環(huán)境搞的軟七八糟的，所以我不想這樣做，所以我決定先建立一個文件夾，來儲存這些輸出文件。所以我該用了下面的語句：#mkdir refrence_index#bowtie2-build Osativa.fa refrence_index/indexname查看一下refrence_index文件下面都有些什么文件。我們會在refrence_index這個目錄下面看到很多個文件，主要為兩種類型，一種是bt2，還有一種是rev.x.bt2。但

11、這些文件都有一個共同特點，文件名的前綴都是indexname，這時你應該能理解bowtie2-bulid這句命令的作用了吧。4.2.利用tophat命令將分別將代比較的reads maping到參考基因組序列上建立完參考基因組序列的索引文件之后，我們就該講樣品的數(shù)據(jù)reads拿去maping了。也許你會問為什么要建立索引？實際上這是為了加快mapping的速度。tophat最簡單的格式是:#tophat -max-intron-length 20000 索引文件 樣品R1文件 樣品R2文件結合本實驗環(huán)境我們可以寫這樣的命令:#tophat -max-intron-length 20000 re

12、frence_index/indexname F1_sample_R1.fastq F1_sample_R2.fastq(備注命令中的“”是用來實現(xiàn)連接下一行的，第一二行本來應該為一個句子)還是那句話，我很討厭輸出的一大堆文件弄的我的實驗環(huán)境亂七八糟的，所以我還是決定為這些輸出文件找個好歸宿，所以我會用下面這幾命令：#tophat -o F1 -max-intron-length 20000 refrence_index/indexname F1_sample_R1.fastq F1_sample_R2.fastq（備注：F1一般就令其為樣品名）考慮到我們還有參考基因組的注釋文件，所以我們可以

13、一起把這個文件加進去。$tophat -o F1 -p 8 -G Osativa.gtf -max-intron-length 20000 refrence_index/indexname F1_sample_R1.fastq F1_sample_R2.fastq可見這句命令使用到了三類文件：l 參考基因組序列的索引文件l 參考基因組的注釋文件l 樣品的數(shù)據(jù)文件，包括R1和R2兩個文件。這里用到了tophat的三個參數(shù)：l -o 指定輸出文件的存放目錄。l -G 指定參考基因組的注釋文件（可有也可沒有，有的話，最后對于分析結果有幫助）。l -p 指定運行該命令使用多少個cpu ,沒有指定情況下

14、默認為一個。 -p 8 則指定使用8個。l -max-intron-length 最大的intron長度。Tophat會忽略長度大于該值的donor/acceptor pairs，除非有l(wèi)ong read支持。根據(jù)實際情況來制定。最后得到的數(shù)據(jù)存放在了F1里面，有不少類型的文件，但最主要的是accepted_hits.bam文件。對F1樣品進行maping操作之后，同樣需要對P1樣本進行同樣的操作，輸出文件制定為P1。$tophat -o P1 -p 8 -G Osativa.gtf -max-intron-length 20000 refrence_index/indexname P1_sa

15、mple_R1.fastq P1_sample_R2.fastq4.3利用cufflinks軟件分別將待測樣品的轉(zhuǎn)錄組reads拼接起來，并同時計算每個樣品各個基因的rpkm值由cufflinks的命令格式Usage: cufflinks options <hits.sam> 可知，該命令只有一類輸入文件，雖然文件類型要求是hits.sam，實際上數(shù)據(jù)hits.bam文件也是可以的，這兩個文件的關系，可通過samtools命令來實現(xiàn)轉(zhuǎn)換。(對于cufflinks不想講太多的東西，因為這個命令輸入文件即我們上一步由tophat產(chǎn)生的F1或P1里的accepted_hits.bam。所

16、以對于樣品F1的操作命令可以寫成：$cufflinks -p 8 -o F1_assemble_transcripts F1/accepted_hits.bamcufflinks參說明：l -o F1_assemble_transcripts 指定輸出文件的存放目錄為 F1_assemble_transcripts 。l -p 8 則指定使用8個cpu跑這個命令。你可以看看F1_assemble_transcripts下都有哪些文件，順便看看最重要的這個genes.fpkm_tracking文件里面的內(nèi)容。自己把下面這個圖片拉大，查看文件內(nèi)容。對樣品P2進行同樣的操作$cd /share/Pu

17、blic/BioinfoTraim/Train_RNAseq$cufflinks -p 20 -o P1_assemble_transcripts P1/accepted_hits.bam4.4.利用coffmerge和cuffdiff軟件計算每個樣品各個基因的fpkm值。這一個步驟中有一些地方會讓一些初學者感到奇怪，比如1. cuffmerge /manual.htmlUsage: cuffmerge Options <assembly_GTF_list.txt>這個命令中的assembly_GTF_list.txt 文件到

18、底指的是什么?實際上，命令幫助里是有說到的“cuffmerge takes two or more Cufflinks GTF files and merges them into a single unified transcript catalog. Optionally, you can provide the script with a reference GTF, and the script will use it to attach gene names and other metadata to the merged catalog”這句話的意思是cuffmerge的作用是將兩

19、個或者更多cufflinks產(chǎn)生的trainscripts.gtf文件合成一個統(tǒng)一的轉(zhuǎn)錄組目錄。既然涉及到兩個或者更多的將兩個或者更多cufflinks產(chǎn)生的trainscripts.gtf文件，那么怎么告訴cuffmerge這個軟件，我這些文件都在哪里呢？我們可以將這些文件所在的地方都存放在一個文本文件里，如慣例使用assembly_GTF_list.txt這個文件來存儲這些信息。我們可以用下面這個語句，來填寫assembly_GTF_list.txt這個文件。$ find -name transcripts.gtf > assemblly_GTF_list.txt2. assembl

20、ly_GTF_list.txt文件準備好之后，我們就開始運行cuffmerge命令吧。$cuffmerge -g Osativa.gtf -s Osativa.fa -p 8 assemblly_GTF_list.txtCuffmerge參數(shù)說明l -o 指定輸出文件存放的目錄，如果你沒指定，它會默認存到merged_asm這個文件里。l -g 可選項，主要是指定參考基因組的注釋文件。 l -s 指定參考基因組的序列文件。現(xiàn)在讓我們來看一下，運行完之后的變化吧，多產(chǎn)生了一個merged_asm目錄，這個文件里有l(wèi)ogs目錄和一個merged.gtf文件。這時候讓我給實驗來上最后沉重的一擊吧,在

21、/share/Public/BioinfoTraim/Train_RNAseq下運行：$ cuffdiff -o F1-vs-P1 -b Osativa.fa -p 20 -L F1,P1 -u merged_asm/merged.gtfF1/accepted_hits.bam P1/accepted_hits.bamUsage:cuffdiff options <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> . sampleN_hits.sam由該命令即可知，cuffdiff需要一下幾個文件l

22、 cuffmerge階段產(chǎn)生的merged.gtf文件l tophat產(chǎn)生的 F1和P1的xx_hits.bam文件cuffdiff參數(shù)說明l -o 指定輸出文件存放的目錄l -b 使用偏定校正，后面接參考基因組的序列，實際作用還不是很清楚。l -p 指定運行該命令使用多少個cpul -L 條件標簽逗號分隔的列表l -u use 'rescue method' for multi-reads 五利用R軟件查看結果文件。1.下載安裝R軟件 2.啟動R軟件，安裝需要的軟件包CummeRbund package.通過輸入以下命令就可以了> source('/biocLite.R')> biocLite('cummeRbund')Use WinSCP to copy folder F1-vs-P1 (/export/xxxx/rnaseq/cufflinks/F1-vs-P1) to your computer. Download it to C: