1  

上圖來源於(http://sourceforge.net/projects/trinityrnaseq/files/misc/RNASEQ_WORKSHOP/rnaseq_workshop_slides.pdf/download)

此套分析流程的開發團隊將此protocol取名為燕尾服(Tuxedo Suite),目的是將NGS技術定序RNA-Seq的reads透過alignment到已知基因體序列後,進行重組轉錄本(transcript reconstruction)、genes/transcripts 定量及分析不同condition間有顯著差異的genes/ transcripts,最後在搭配R語言將統計結果圖像化;此套分析流程適用於有已知基因體序列的物種(例如:人類、老鼠、斑馬魚、阿拉伯芥、番茄..等)。

  2  

上圖來源於( http://sourceforge.net/projects/trinityrnaseq/files/misc/RNASEQ_WORKSHOP/rnaseq_workshop_slides.pdf/download )

 

燕尾服的開發團隊作者如下圖,是一個橫跨多所大學及研究機構所完成的分析protocol。

3  

上圖來源於(http://sourceforge.net/projects/trinityrnaseq/files/misc/RNASEQ_WORKSHOP/rnaseq_workshop_slides.pdf/download)

 

由於RNA-Seq的reads會有包含橫跨exon-junction的情形,燕尾服protocol的第一個步驟是使用tophat將reads map回genome,如下圖,概念上簡單說是當reads無法end-to-end的map回genome時,這種reads會被切為多個片段(segment),透過小片段的segment map回genome後搭配flanking找出跨越exon的reads。

4  

上圖來源於(http://genomebiology.com/2013/14/4/R36/figure/F6)

透過tophat將reads map回genome後(下圖a),接著會透過Cufflinks來重建轉錄本(Transcript Reconstruction),這個步驟是以genome為參考序列進行assembly (reference-based assembly);由於pre-mRNA會透過alternative splicing形成不同的transcript isoforms, 透過橫跨exon-junction的reads能回推有哪些類型的transcript isoforms存在(下圖b,c,d)。

5  

6  

上圖來源於(http://sourceforge.net/projects/trinityrnaseq/files/misc/RNASEQ_WORKSHOP/rnaseq_workshop_slides.pdf/download)

 

經過transcript reconstruction後,將進行基因及轉錄本的定量(quantification),定量的概念如下圖,當我們將reads map到genome後,可得知每一個gene/transcript有多少條 reads覆蓋 (稱為read count),可以想像越長的基因,能被覆蓋到的read count就越多,所以會將read count normalize 基因的長度及本次實驗的定序量,定量的單位有RPKMFPKM,在Cufflinks的方法中,是使用FPKM為基因表現量的單位。

7  

上圖來源於(http://sourceforge.net/projects/trinityrnaseq/files/misc/RNASEQ_WORKSHOP/rnaseq_workshop_slides.pdf/download)

 

當定量完後,若實驗設計有多個組別要互相比較表現量,則可透過Cuffdiff計算顯著表現的genes/ transcripts及使用CummeRbund呈現統計結果。

 

簡單介紹燕尾服Protocol的基本概念後,接著要教學如何操作這個protocol。本次要練習的資料,實驗設計是假設有兩個不同Condition的RNA-Seq資料,我們想找出這兩個condition下最後有顯著差異的基因有哪些,且本次實驗是有已知的基因體參考序列。流程圖如下圖

8   

上圖來源於(http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html)

 

步驟一: 安裝 VirtualBox

請先安裝VirtualBox,用來建立虛擬電腦,安裝作業系統

Windows載點

OS X載點

 

這邊以Window 7 下的執行當作範例,所有的步驟都按照預設即可。

9  

 

步驟二: 下載虛擬電腦的image 檔並匯入VirtualBox

之後請下載Trinity 團隊提供的Ubuntu image

此虛擬電腦為linux Ubuntu的作業系統,已經內建許多NGS常用的軟體。

包含BowtieBowtie2Tophat2CufflinksCummeRbundSamtoolsigvTrinitygmapR…等。

對於初次接觸生物資訊的初學者來說,可能連linux的作業系統都沒用過,但生物資訊多數開發的軟體是以linux的平台為基礎,此練習省略許多安裝上的門檻,可以直接練習分析RNA-Seq資料常用到的軟體及技巧。

 

VirtualBox安裝完成後請打開

 10  

 

選擇檔案->匯入應用裝置->選擇剛剛下載的Ubuntu_RNASeqWorkshop.ova檔。

11  

 

12  

 

13  

 

等待匯入。

14  

 

當匯入完成後,可看到虛擬電腦已經建立完成,可以看到預設的基礎記憶體是4096MB,若你們練習的電腦RAM沒有那麼大,請點”設定值”,將此值調到1024MB,避免開起虛擬電腦時整個電腦Lag。

15  

 

接者請案”啟動”,則會開啟虛擬電腦。

16  

 

如點選啟動後出現以上的錯誤訊息,請重開機進入BIOS,將Intel Virtualization Technology 開啟,下圖是以筆者的主機板BIOS為例。

17  

 

步驟三: 開啟虛擬電腦,並打開Terminal

開機完成後的畫面如下,此為Ubuntu的桌面。

18  

 

請從左側的選項裡找到Terminal並開啟。

 

開啟終端機的畫面如下。

20  

 

步驟四: 進入rnaseq_workshop_data的資料夾並解壓縮範例的壓縮檔

在開始練習前先簡單說明使用 terminial 時,使用 pwd 指令可以看到目前的資料夾路徑 ( 開啟 terminal 的初始路徑為 /home/ubuntu) ,要切換目錄會使用到 cd 這個指令,若要看目前這個目錄有甚麼檔案,則使用 ls -alh ,在畫面操作的結果如下圖。

> cd rnaseq_workshop_data

> ls -alh

21  

  如上圖可以看到rnaseq_workshop_data裡有一些範例檔,簡單說明有condition A的Paired-end reads(condA.left.fa.gz及condA.right.fa.gz)、condition B的Paired-end reads(condB.left.fa.gz及condB.right.fa.gz)、物種的參考序列(genome.fa.gz及bowtie和bowtie2的index檔案)

 

因為這些範例檔有些有先被壓縮,副檔名為gz,在此使用gunzip這個指令來解壓縮rnaseq_workshop_data這個資料夾內的壓縮檔,如下圖執行完後請在使用ls指令可看到副檔名為gz的檔案皆已經被解壓縮。

> gunzip *.gz

22  

 

步驟五: 針對Condition A 及B 執行Tophat 、Samtools 及Cufflinks

接著請依序執行下列三個指令,針對condition A的reads執行tophat及Cufflinks,執行的情形請見下圖;-l參數為最大的intron長度 -i為最小的intron長度 -o的參數為產生的report的資料夾,result的檔案細節及格式請參見(http://tophat.cbcb.umd.edu/manual.shtml)。

> tophat –I 1000 –i 20 –o condA_tophat_out genome condA.left.fa condA.right.fa

> samtools index condA_tophat_out/accepted_hits.bam

> cufflinks -o condA_cufflinks_out condA_tophat_out/accepted_hits.bam

   23  
     

針對condition B也依此類推,請依序執行下面三行指令,執行的畫面請見下圖

> tophat –I 1000 –i 20 –o condB_tophat_out genome condB.left.fa condB.right.fa

> samtools index condB_tophat_out/accepted_hits.bam

> cufflinks -o condB_cufflinks_out condB_tophat_out/accepted_hits.bam

 24  

 

 

步驟六: 執行cuffmerge

以下4個指令是將condition A及 B 使用cuffmerge將結果merged起來,執行後會有merged_asm的資料夾產生,merged.gtf為結果,此4個指令的執行情形請見下圖。

> echo condA_cufflinks_out/transcripts.gtf > assemblies.txt

> echo condB_cufflinks_out/transcripts.gtf >> assemblies.txt

> cat assemblies.txt

> cuffmerge -s genome.fa assemblies.txt

25  

 

步驟七: 使用IGV 瀏覽alignment 資訊

完成上述的步驟時,此時的檔案路徑在home/ubuntu/rnaseq_workshop_data,請先回到/home/ubuntu,並且使用unzip解壓縮IGV的壓縮檔,並進入IGV_2.3.8的資料夾內開啟IGV程式。

> cd ../

> unzip IGV_2.3.8.zip

> cd IGV_2.3.8

> sh igv.sh

26  

 

開啟igv後,請選Genomes-> Load Genome from File -> 選擇genome.fa

27  

 

接著請選File-> Load from File -> 依序選擇condA_tophat_out/accepted_hits.bam、condB_tophat_out/accepted_hits.bam、merged_asm/merged.gtf、genes.bed

28  

 

完成後可觀看這次conA及conB的alignment情形,如下圖,瀏覽IGV的方法請參考(http://www.broadinstitute.org/igv/Navigate)。

29  

 

步驟八: 執行cuffdiff 挑出有顯著差異的genes/transcripts

接著請關掉IGV回到Terminal,此時的路徑為/home/ubuntu/IGV_2.3.8,請再次回到/home/ubuntu/rnaseq_workshop_data/,並執行cuffdiff,結束後的結果在diff_out資料夾下,result的檔案內容細節請參見(http://cufflinks.cbcb.umd.edu/manual.html#cuffdiff_output),執行的情形如下圖。

> cd ../rnaseq_workshop_data

> cuffdiff -o diff_out -b genome.fa -L condA,condB -u merged_asm/merged.gtf condA_tophat_out/accepted_hits.bam condB_tophat_out/accepted_hits.bam

30  

 

步驟九: 使用R 語言匯入CummeRbund library 及Cuffdiff 的結果

R語言是一個統計上常使用的程式語言,在linux中要開啟最簡單的作法是在終端機打”R”,即可進入R的編輯環境,如下圖。

> R

31  

 

接著請匯入CummeRbund的library,執行的情形如下圖。

> library(cummeRbund)

32  

 

再來是將步驟8執行Cuffdiff的結果匯入到R中,執行的情形如下圖

> cuff=readCufflinks(‘diff_out’)

33  

 

步驟十: 使用CummeRbund 產生統計圖表

CummeRbund能以Density plotScatter plotVolcano plot..等 圖表來呈現多個sample內的基因表現量情形,其指令及產生的圖表如下。

 

Examine the distribution of expression values for the reconstructed transcripts

> csDensity(genes(cuff))

 

Examine transcript expression values in a scatter plot

> csScatter(genes(cuff),'condA','condB')

 

Examine individual densities and pairwise scatterplots together.

> csScatterMatrix(genes(cuff))

 

Volcano plots are useful for identifying genes most significantly differentially expressed

> csVolcanoMatrix(genes(cuff),'condA','condB')

 

Extract the‘genes’that are significantly differentially expressed

下列第二航的指令執行後,可看出有顯著差異的基因的筆數,如下圖,共有436筆,第三航指令是將p-value < 0.1的挑出來,可看出有74筆

> gene_diff_data=diffData(genes(cuff))

> nrow(gene_diff_data)

> sig_gene_data = subset(gene_diff_data,(significant=='yes' | p_value < 0.1))

38  

 

此步驟的基因是以使用者挑選自行興趣的基因,在此以XLOC_000006為例子,可見上圖XLOC_000006為P-value < 0.1的基因。

> var_XLOC_000006=getGenes(cuff,’XLOC_000006’)

39  

 

> expressionBarplot(var_XLOC_000006,logMode=T,showErrorbars=F)

Plot the gene expression values for the gene under each condition

 

參考資料:

  1. Trapnell C,Adam Roberts .,et al.   Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature Protocols 7,562–578 (2012) doi:10.1038/nprot.2012.016
  2. Trapnell C, Pachter L .,et al. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics doi:10.1093/bioinformatics/btp120
  3. Kim D, Pertea G .,et al. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. . Genome Biology 2011, 14:R36
  4. Trapnell C, Williams BA .,et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology doi:10.1038/nbt.1621
  5. Roberts A, Trapnell C .,et al. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biology doi:10.1186/gb-2011-12-3-r22
  6. Roberts A, Pimentel H .,et al. Identification of novel transcripts in annotated genomes using RNA-Seq. Bioinformatics doi:10.1093/bioinformatics/btr355
  7. Trapnell C, Hendrickson D .,et al. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nature Biotechnology doi:10.1038/nbt.2450

 

 

 

 

logo yourgene  

YourGene 發表在 痞客邦 PIXNET 留言(0) 人氣()