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

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

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

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

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

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

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

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

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

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

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