2. 程式人生 > >使用Tophat+cufflinks分析差異表達

使用Tophat+cufflinks分析差異表達

阿新 • • 發佈：2018-11-09
使用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 做一些常見的統計圖

 

相關推薦

使用Tophat+cufflinks分析差異表達

使用Tophat+cufflinks分析差異表達  2017-06-15 19:09:43     522     0     0   &nbs

使用Trinity拼接以及分析差異表達一個小例子

使用Trinity拼接以及分析差異表達一個小例子  2017-06-12 09:42:47     293     0     0 Trinity 將測序資料

轉錄組差異表達分析小實戰（一）

轉錄組差異表達分析小實戰（一）  Posted: 七月 28, 2017  Under: Transcriptomics  By Kai  no Comments 讀文獻獲取資料 文獻名稱：AKAP95 r

轉錄組差異表達分析小實戰（二）

轉錄組差異表達分析小實戰（二）  Posted: 八月 14, 2017  Under: Transcriptomics  By Kai  no Comments 差異基因表達分析 我按照前面的流程轉錄組差

使用GEO資料庫來篩選差異表達基因，KOBAS進行KEGG註釋分析

前言 本文主要演示GEO資料庫的一些工具，使用的資料是2015年在Nature Communications上發表的文章Regulation of autophagy and the ubiquitin-proteasome system by the FoxO transcriptional networ

差異表達基因分析[轉載]

轉自：https://wenku.baidu.com/view/2532ab5176c66137ef06191a.html  1.轉錄組 2.轉錄組研究重要性 3.轉錄組研究技術 3.1三種的比較 3.2技術比較 4.高通量測序技術的應用 5.轉錄組

利用R語言對RNA-Seq進行探索分析與差異表達分析

介紹 本文參考 bioconductor 中RNA-Seq workflow: gene-level exploratory analysis and differential expression並對其根據需要進行了增減。 試驗資料 資料

生物資訊學入門 使用 GEO基因晶片資料進行差異表達分析（DEG）——Limma 演算法 資料 程式碼 結果解讀

       差異表達分析通常作為根據基因表達矩陣進行生物資訊學分析的第一步，有助於我們觀察基因在不同樣本中的表達差異，從而確定要研究的基因和表型之間的聯絡。常用的基因表達資料來自基因晶片或高通量測序。雖然矩陣看起來差不多，但是由於服從不同的分佈，因此在進行差異表達的時候需要

單細胞分析實錄(7): 差異表達分析/細胞型別註釋

> 前面已經講解了： > [單細胞分析實錄(1): 認識Cell Hashing](https://www.jianshu.com/p/02b70bda5d0a) > [單細胞分析實錄(2): 使用Cell Ranger得到表達矩陣](https://www.jianshu.com/p/5a

RNA-seq中的基因表達量計算和表達差異分析

RNA-seq中的基因表達量計算和表達差異分析 差異分析的步驟：1）比對；2） read count計算；3） read count的歸一化；4）差異表達分析； 背景知識：1）比對：普通比對： BWA，SOAP開大GAP比對：Tophat（Bowtie2）；2） Read count(多重比對的問題）：丟

常用正則表達式爬取網頁信息及HTML分析總結

logfile mpi 開始 order 標題 ear 爬取網頁 常用 enter Python爬取網頁信息時，經常使用的正則表達式及方法。 1.獲取<tr></tr>標簽之間內容 2.獲取<a href..></a>超鏈接

Python日誌分析與正則表達式

logs sea 篩選 ear d+ class 時間 針對 日誌 程序員經常會面臨日誌的分析工作。而正則表達式是處理日誌的必備工具。 “Line 622: 01-01 09:04:16.727 <6> [pid:14399, cpu1 dabc_pwym_t

Python基礎知識進階（五---2）----程序基本結構、簡單分支、異常處理、三大實例分析、基本循環結構、通用循環構造方法、死循環嵌套循環、布爾表達式

方法 算法 嵌套 構造方法 決策樹 輸入 繼續 實例 控制 上一篇隨筆寫的內容有點多了，決定分成兩節，不然自己看的時候也頭疼。 三者最大實例: 分支結構可以改變程序的控制流，算法不再是單調的一步步順序執行。 　　　　假設：以找出三個數字中最大者的程序設計為例。 　　　　　

逗號表達式的分析

直接 計算 表達 返回 逗號運算符 每次 什麽 賦值 color fun函數返回的a值和傳入的a值無關，每次循環得到的k值和上一次循環也無關，所以直接看最後一次循環結果就好。 結果的關鍵在於逗號表達式 a=(c++,b++); 看了大家的解析，想要提一點，以下兩個表達式結果

日誌分析代碼實現(正則表達式)

日誌日誌分析代碼實現(正則表達式)思路 使用正則表達式處理: 使用正則提取對應內容 每段數據轉換為對應格式 精簡代碼,異常處理,代碼效率檢查import datetimeimport relogline = ‘‘‘183.60.212.153 - - [1

C#正則表達式的遞歸匹配分析

重復字符 renren href ffffff share lock src www. net 在C#程序設計中經常會遇到這樣的需求，要求匹配出成對的小括號裏的內容，但是一般正則表達式中的 ?R 的語法似乎在C#中不被支持， 經過一番查找與測試，終於找到以下一段描述 /

Python數據分析學習-re正則表達式模塊

如果 dal nds 隨機數 choice 而且 inf groups 對象 正則表達式 為高級的文本模式匹配、抽取、與/或文本形式的搜索和替換功能提供了基礎。簡單地說，正則表達式（簡稱為 regex）是一些由字符和特殊符號組成的字符串，它們描述了模式的重復或者表述多個字

不同系統裏同一Customizing activity的顯示差異分析

custom sys sin 查看 dmi inf 二維 開關 admin 比如SAP HANA Live Reporting這個Customizing Activity，在System AG3裏不可見，但是在另一個系統QHD裏卻能使用。比較下列兩張圖的差異。 一種可能

5個適用於初學者的有用數據分析表達式（DAX）函數

src 總數 圖片 表達 scene 創建 之間 top n graph 數據分析表達式（DAX）入門可能令人生畏，但是，當你了解一些基本功能後，就可以幫助你解答許多有關數據的新見解。雖然在Power BI或Pivot Charts中創建視覺效果很容易，但我們經常希望

原方法與正則表達式方法對比性能分析

pre 帶來 mil .com 數字 結對編程 前言 興趣 誤差 0．前言 　　寫這篇隨筆的原因呢，還得感謝此次結對編程與南通大學鞠老師的合作，因時間較短開發出來的小程序在性能上並沒有足夠多的考慮。其實我一直以來都有一種習慣，就是對於一些細枝末節的程序，並不太願意花時間去考