碼上快樂

相關內容    簡體    繁體

使用Tophat+cufflinks分析差異表達

本文轉載自   查看原文   2018-11-09 22:17   2756    生物信息學/ 轉錄組Transcriptome/ RNA-seq
使用Tophat+cufflinks分析差異表達
 2017-06-15 19:09:43     522     0     0

 

                           

 

使用TopHat+Cufflinks的流程圖

 

 

 

    序列的比對是RNA分析流程中核心的一步。序列的比對，或者說是字符串的比對本身就是計算機科學中的一個經典問題，在生物信息學中更加頻繁的出現。序列比對中的錯配，插入、缺失可以識別出樣本和基因組之間的多態性，甚至可以找出腫瘤樣本中的gene fusion。而map到沒有注釋的基因可能是新的編碼基因，或者是非編碼RNA。同時RNA-seq的序列比對可以揭示新的選擇性剪接和同工型（isoform）。

    此外，序列的比對也可以用作精確定量基因或者轉錄本的表達量，因為顯然表達豐度與產生的reads是直接成比例的（需要標准化）

    Tophat使用Bowtie作為比對引擎，Bowtie使用了一種極其緊湊的數據結構FM index 來儲存參考基因組的序列，並且能夠迅速的查找（ tens of millions per CPU hour）。但是Bowtie的比對不允許gap的出現（Bowtie2中已經可以允許了），所以Bowtie不能比對跨越了內含子的reads。

    Tophat將Bowtie不能align的reads打斷成更小的片段稱作segments。一般情況下，當單獨處理時，這些segmens可以map到基因組，當一個read的幾個segmens可以map到基因組時,Tophat推斷這個read跨越了剪接位點同時推測剪接位點的位置。通過處理每個‘initially unmappable’ read，Tophat在沒有剪接位點注釋的信息下能夠在轉錄組上建立起一個剪接位點的索引。但是這個剪接位點圖譜的構建還不夠完整，即使是在一些深入研究的模式植物中，每個轉錄組測序中都能發現新的剪接位點。

Transcript assembly with Cufflinks

    計算基因表達量的另一個問題是，因為選擇性剪接的原因，幾個不同的轉錄本（isoform）可能擁有相同的外顯子，此時難以確定reads到底來自其中哪個轉錄本（isoform）。所以能否確定所有的splice variants（isoform）決定着表達量計算的准確性。而這個又很難確定，所以cufflinks通過map到基因組的reads組裝起一個簡陋的轉錄組，用reads拼接成含有重疊部分但是長度不同的轉錄本（稱作”transfrags“，作為splice variants的推測。拼接以后，Cufflinks計算使用嚴格的統計模型來計算每個transfrag的表達量。

    當有多個樣本的時候，一種方法是將所有樣本的reads合起來，拼接成一個轉錄組。這種方法的缺點是：

  1. 大量reads帶來的計算不便，需要更高配置的服務器和大量的時間（可以參考 de no vo Trinity，動輒上百G的內存需求）
  2. 多個樣本重疊使得確定所有的splice variants（isoform）更加困難

所以Cufflinks采用的策略是先單獨拼接每一個樣本的reads，然后使用cuffmerge來綜合所有樣本的拼接結果

    

 Differential analysis with Cuffdiff

Cufflinks還包括一個cuffdiff程序，用於計算基因表達豐度和統計推斷。

除了基本的差異表達分析外，cuffdiff還

 

 

運行環境需求：

  • Bowtie2
  • SAM tools         

要求Bowtie2和SAM tools 都已經安裝且添加到系統環境變量中

  • Tophat
  • Cufflinks
  • Cummerbund
  • 64-bit linux or Mac os x，最低要求是4G內存，推薦16G內存以上

 

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0125031#references

數據來自於plos one上的一篇文獻，測序策略是一小部分測序數據用來拼接轉錄組（雙端），因為當時茶樹的基因組還沒有公布，只能deno vo拼接。剩下的做定量的單端測序，兩個重復。

今年發表了一篇茶樹基因組的文章，上傳了一個拼接的基因組fasta文件，我們可以用做參考。

TopHat

安裝

Tophat在運行中會調用Bowtie1 或者Bowtie2，所以首先要確定你的系統中安裝了Bowtie，並且已經添加到了環境變量中

Tophat提供了已經編譯好的包可以直接下載，也可以直接下載源碼編譯，為了方便我們當然選擇前一種。網頁下載好壓縮包或者直接

  1. wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.1.1.Linux_x86_64.tar.gz
  2. tar -xvf tophat-2.1.1.Linux_x86_64.tar.gz
  3. export PATH="${PATH}:/usr/local/biosoft/tophat-2.1.1.Linux_x86_64/"
  4. #添加快捷方式
  5. ln -s ~/tophat-2.0.0.Linux_x86_64/tophat2 .

解壓后將目錄添加到環境變量中，或者添加一個快捷鏈接

為了確保安裝正確，最好進行一下測試，尤其是使用編譯源碼的方式安裝的時候。下載test_data

  1.     
  2.     tar -xvf test_data.tar.gz
  3.     cd test_data
  4.     tophat -r 20 test_ref reads_1.fq reads_2.fq

結果也沒讓我們失望，果然失敗了。。

研究了一下原因，發現是因為bowtie2的版本太新的原因，換成比較舊的版本就ok了

tophat 官網上最新更新的版本還是2016/02/23的，而且官網上也強調了

Please note that TopHat has entered a low maintenance, low support stage as it is now largely superseded by HISAT2 which provides the same core functionality (i.e. spliced alignment of RNA-Seq reads), in a more accurate and much more efficient way.

 建立基因組索引

對於一些研究比較深入的物種，可以下載已經建立好的索引 using a pre-built index

或者根據你自己的fasta參考文件使用Bowtie2建立，

  1. bowtie2-build /home/shady/RNA-seq/reference_genome_tea/Teatree_Assembly.fas Tea_genome

運行成功后產生后綴為.1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2, and .rev.2.bt2  這些文件組成了索引

 

用我自己的筆記本，大概3G的基因組，build了將近一個小時。。

 

准備reads

Tophat 可以接受FASTQ or FASTA 格式，推薦fastq

可以用fastqc看一下數據的質量，主要是看看有沒有特別差的數據以及接頭有沒有去除

fastqc有圖形界面，也可以使用命令行模式

具體的參數和結果輸出請看這篇文章，講的比較詳細

因為是已經發表了的數據，這里只是簡單的看一下

  1. #直接安裝
  2. sudo apt-get fastqc
  3. fastqc *.gz --outdir=/home/shady/RNA-seq/quality_control --extrct
  4. #結果對每個文件都生成一個html文件和一個壓縮包 
  5. #參數--extrcat 會自動解壓壓縮包，進去看一下summary.txt
  6. cd /home/shady/RNA-seq/quality_control/SRR1747100-25du-4h-1.sra_fastqc
  7. cat summary.txt

運行 Tophat 

Tophat首先會調用Bowtie2使用全局比對的方法map你的reads到參考基因組，因為你的reads來自於經過剪接的轉錄本，所以reads  map到基因組位置堆積起來就像一個參考基因組上的"islands"，許多個這樣的"island" 就組成了所謂的外顯子

Tophat然后運行一個程序，使用沒有map到的reads來尋找剪接位點

一個典型的使用方法為

使用如下腳本處理我們的6個數據

  1. #參數 -p 表示運行的線程數
  2. #參數 -G 可以選擇基因組的注釋文件（可選）
  3. ls `pwd` |grep "gz"| while read filename
  4. do 
  5. tophat -p 4 -G Teatree.gff3 -o ${filename%%.*} Tea_genome ${filename}
  6. done

參數詳解  一般用默認參數

輸出

Tophat運行結束后，將在當前目錄產生以下幾個文件，可以直接在基因組瀏覽器中查看如IGV等

  1. accepted_hits.bam   一個SAM格式的read alignment 列表，SAM是一種緊湊的短序列比對格式文件，里面包含可以map到基因組上的比對信息
  2. junctions.bed       一個UCSC BED 文件。是ucsc 的genome browser的一個格式，報告剪接位點 
  3. insertions.bed and deletions.bed 

 

 

參數解釋

 常用

  1. o 輸出目錄，默認值為 “./tophat_out”。  

  2. –solexa-quals/solexa1.3-quals 質量編碼，關於質量編碼格式請參考《FastQ格式介紹》  

  3. -p 線程數，默認值為單線程1.，可以使用多線程  

  4. -G/–GTFSupply TopHat with a set of gene model annotations and/or known transcripts, as a GTF 2.2 or GFF3 formatted file.指定已有轉錄本信息  

  5. –no-novel-juncs 不查找新的可變剪切  

  6. -r 比對時兩成對引物間的距離中值。比如說，如果你的插入片段有300bp，而每個引物有50bp，那么r值就應該是200=(300+50*2)/2。沒有默認值，如果是末端配對比對時這個值是必須的。  

  7. –mate-std-dev 末端配對時中間插入片段的長度的標准差，默認值為20bp  

 參考：http://ccb.jhu.edu/software/tophat/manual.shtml#output

 

使用Cufflinks計算差異表達

官方手冊 

Cufflinks包含了以下程序：

  • Cufflinks 使用經過Tophat比對好的SAM格式來拼接轉錄組（每一個樣本都進行拼接），如果有基因組注釋文件的話，效果會更好。這個步驟的目的在前面已經講過，包括后面的cuffmerge，即為了更精確地確定轉錄組的所有轉錄本以及其所有的剪接變體，從而更精確的進行定量分析。

        如果已經有了參考基因組和基因組注釋文件，這個步驟和cuffmerge都可以跳過，也就是說      

        cuffquant/ cuffdiff 可以直接接受Tophat輸出的文件進行定量分析，差別就是選項中的基因

        組注釋文件用你之前下載的基因組注釋文件，而不是由cufflinks、cuffmerge生成的注釋文

        件，一般不推薦這么做。因為即使是研究比較深入的模式植物，每次轉錄組測序都會發現新

        的轉錄本。

  • cuffmerge 將Cufflinks拼接的所有轉錄組整合起來
  • cuffquant v2.2.0 后出現的新功能，在之前的cuffmerge 和 cuffdiff 之間插入了一步，作用是計算每個樣本的表達量（應該是raw count），產生一個中間文件，這個文件可以用cuffdiff  和cuffnorm分析，是一個可選步驟。
  • cuffdiff 可接受Tophat的輸出/cuffquant的輸出/或者跳過cuffquant

         對每個處理都兩兩進行差異分析，分析哪些基因上調或者下調，顯著性如何等，后面介紹。

  • cuffnorm 也是v2.2.0 后新出來的功能，不進行差異分析，只計算FPKM等值，進行標准化

 

下面腳本進行了從Cufflinks到cuffdiff的一個流程，絕大部分參數使用默認參數即可，其他參數請看官方手冊，因為中間涉及幾個不同程序，需要要搞清楚每個程序需要的輸入文件和輸出文件。

 

  1. #運行Cufflinks，分別拼接每一個樣本的轉錄組
  2. ls `pwd` |grep 'gz'|while read filename
  3. do 
  4. cufflinks -p 4 -o "./${filename%%.*}_clout" "${filename%%.*}/accepted_hits.bam"
  5. done
  6. #創建一個 assemblies.txt
  7. #運行cuffmerge，將所有樣本拼接好的轉錄組合成一個
  8. #參數 -g 添加基因組注釋文件，可選
  9. #參數 -s 為參考基因組 fasta格式
  10. cuffmerge -g Teatree.gff3 -s Teatree_Assembly.fas -o ./cuffmerge_out/-p 4 assemblies.txt
  11. #運行 cuffquant
  12.  
  13. ls `pwd` |grep 'gz'|while read filename
  14. do
  15. cuffquant -p 4 -o "./cuffquant_result/${filename%%.*}_cqout/" ./cuffmerge_out/merged.gtf  "${filename%%.*}/accepted_hits.bam"
  16. done
  17.  
  18. #運行cuffdiff
  19. cuffdiff -p 4 -o ./cuffdiff_out -L T25_4h,T4_4h,T4_8h ./cuffmerge_out/merged.gtf\
  20.  ./cuffquant_result/SRR1747100-25du-4h-1_cqout/abundances.cxb\
  21. ,./cuffquant_result/SRR1747101-25du-4h-2_cqout/abundances.cxb\
  22.  ./cuffquant_result/SRR1747102-4du-4h-1_cqout/abundances.cxb\
  23. ,./cuffquant_result/SRR1747103-4du-4h-2_cqout/abundances.cxb\
  24.  ./cuffquant_result/SRR1747105-4du-8h-1_cqout/abundances.cxb\
  25. ,./cuffquant_result/SRR1747107-4du-8h-2_cqout/abundances.cxb\
  26. ~                                                              

最終cuffdiff運行結束后產生的結果如下

產生了一大堆文件，大致分為以下幾類

本次分析有3個樣本每個樣本兩個重復，共6個文件

FPKM trcking files 計算的每個樣本的FPKM值

其中又分為如下四種，下面的類似

isoforms.fpkm_tracking Transcript FPKMs
genes.fpkm_tracking Gene FPKMs. 因為一個基因可以有多個轉錄本，所以gene FPKM值是具有相同基因id的轉錄本的Fpkm值的總和
cds.fpkm_tracking Coding sequence FPKMs. 具有相同蛋白id轉錄本fpkm值的總和
tss_groups.fpkm_tracking Primary transcript FPKMs. 具有相同轉錄起始位點的轉錄本FPKM值的總和

 

Count tracking files  計算的每個樣本的raw ount數，依舊分為

       

isoforms.count_tracking Transcript counts
genes.count_tracking Gene counts. Tracks the summed counts of transcripts sharing each gene_id
cds.count_tracking Coding sequence counts. Tracks the summed counts of transcripts sharing each p_id, independent of tss_id
tss_groups.count_tracking

Primary transcript counts. Tracks the summed counts of transcripts sharing each tss_id

Read group tracking files  計算的樣本里每個重復的表達數 和count數

isoforms.read_group_tracking Transcript read group tracking
genes.read_group_tracking Gene read group tracking. Tracks the summed expression and counts of transcripts sharing each gene_id in each replicate
cds.read_group_tracking

Coding sequence FPKMs. Tracks the summed expression and counts of transcripts sharing each p_id, independent of tss_id in each replicate
tss_groups.read_group_tracking

 

Primary transcript FPKMs. Tracks the summed expression and counts of transcripts sharing each tss_id in each replicate

Differential expression tests

樣本之間的兩兩比較的差異表達和統計推斷

isoform_exp.diff Transcript-level differential expression.
gene_exp.diff Gene-level differential expression. Tests differences in the summed FPKM of transcripts sharing each gene_id
tss_group_exp.diff Primary transcript differential expression. Tests differences in the summed FPKM of transcripts sharing each tss_id
cds_exp.diff Coding sequence differential expression. Tests differences in the summed FPKM of transcripts sharing each p_id independent of tss_id

header的具體信息

Differential splicing tests - splicing.diff

對每個轉錄前體產生的可變剪接變體的數目

Differential coding output - cds.diff

 在這個文件中只記錄了能產生多個CDS的基因的差異表達數

Differential promoter use - promoters.diff

在這個文件中只記錄了能產生不同轉錄前體（i.e. muti-promoter genes）的基因的差異表達

Read group info - read_groups.info

記錄了樣本重復之間計算定量的一些信息

 

使用cummeRbund 做一些常見的統計圖

 


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 tophat cufflinks cuffcompare cuffmerge 的使用 使用bowtie2 tophat2 及 cufflinks 處理RNA-SEQ數據 轉錄組差異表達分析小實戰（一） 單細胞測序數據的差異表達分析方法總結 R語言——limma包進行多組差異表達分析 簡單使用limma做差異分析 limma package計算差異表達 tophat的用法 鴻蒙輕內核源碼分析：掌握信號量使用差異 Excel-漏斗圖分析（差異分析）
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM