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

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

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

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

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

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

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

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

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

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

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