是時候來一波RNA-Seq差異表達分析實操了

最新 06-16

案例簡介

數據來自於發表在Nature commmunication 上的一篇文章「Temporal dynamics of gene expression and histone marks at the Arabidopsis shoot meristem during flowerin」。原文用RNA-Seq的方式研究在開花階段,芽分生組織在不同時期的基因表達變化。

原文的流程是： TopHat -> SummarizeOverlaps -> Deseq2 -> AmiGO

其中比對的參考基因組為TAIR10 ver.24 ，並且屏蔽了ribosomal RNA

regions (2:3471–9557; 3:14,197,350–14,203,988)。Deseq2隻計算至少在一個時間段的FPKM的count > 1 的基因。

數據存放在http://www.ebi.ac.uk/arrayexpress/, ID為E-MTAB-5130。

實驗設計： 4個時間段（0,1,2,3），分別有4個生物學重複，一共有16個樣品。

數據下載

首先下載數據說明文件：

然後根據數據說明文件提供的FTP鏈接下載

根據下載速度，你們可以選擇去吃吃飯，還是睡睡覺。

下載完RNA-Seq數據後，我們還需要下載一個擬南芥cDNA序列(記住是轉錄組，而不是全基因組，好了，你別說了，我記住了)

然後用Salmon建立索引：

數據定量

由於樣本一共有16個，不可能一條條輸入命令，所以我們寫一個腳本：

根據你電腦的配置，你可以選擇吃下午茶，還是選擇睡個午覺。

數據導入R

當然你完全不必真的去睡午覺，我們可以程序運行的時候準備一下所需要的輸入文件。

可以糾正不同樣本基因長度的潛在改變（比如說differential isoform usage）；能夠用於導入 (Salmon,Sailfish,kallisto)程序的輸出文件；能夠避免丟棄那些比對到多個基因的同源序列，從而提高靈敏度

In particular, thetximportpipeline offers the following benefits: (i) this approach corrects for potential changes in gene length across samples (e.g. from differential isoform usage) (Trapnell et al. 2013), (ii) some of the upstream quantification methods (Salmon,Sailfish,kallisto) are substantially faster and require less memory and disk usage compared to alignment-based methods that require creation and storage of BAM files, and (iii) it is possible to avoid discarding those fragments that can align to multiple genes with homologous sequence, thus increasing sensitivity (Robert and Watson 2015)

雖然的參數看起來很多，但其實需要我們準備的就是兩個和

存放的是salmon的輸出文件，所以我們需要根據文件存放位置，進行聲明

吐槽: 原本的我以為sample從1到16會和194到206是一一對應，但是萬萬沒想到，他的對應關係居然是下圖這個樣子的。看來凡是都不能想當然，一定要仔細看他們的對應關係。

然後我們還要準備一個基因名和轉錄本名字相關的數據框

如果你電腦跑的夠快，基本上這個時候就可以導入數據了。

由於後續要用DESeq2包做差異表達分析，所以需要用這個函數，當然你還需要說明你的實驗設計

數據預過濾：

差異表達分析

標準的差異表達分析步驟在DESeq2隻需要一個函數就可以完成。然後可以通過提取結果表，包括log2 change, p values, adjusted p values.通過查看描述性結果

results還有許多參數：

用於指定比較的對象，默認是最後一個（實驗組）與第一個（對照組）比較。

: 用於指定log2 fold change的閾值，默認是0. fold change = (group1 - group2)/min(group1,group2).

: 顯著性水平閾值

：多重實驗校正

於是我分別用Day1,2,3與Day0進行比較：

最後找到了208個差異表達基因，105上調，103個下調；而文章是298個，148上調，156下調。

探索性分析

特殊基因的表達情況

有些基因是分生組織的特異基因，如STM, KNOTTED-LIKE, FROM ARABIDOPSIS THALIANA 1(KNAT1), CLAVATA 1(CLV1)和CLV3, 從理論上應該有一定的表達量

早花同源異形基因（homeotic gene)，如AP1, APETALA 3(AP3), CAULIFLOWER(CAL)和發育中期和後期基因JAGGED(JAG), WUSCHEL RELATED HOMEOBO 1(WOX1) WOX3從理論上應該不表達。

上面的基因名我們需要先轉成TAIR的ID，才能從dds中提取相應的counts.然後使用和組織數據, 以meristem_identify_genes為例：

由於這個畫圖過程，後續會用大很多次，為了避免重複操作，我寫成了一個函數，我建議你根據上面的流程自己寫一個函數，不要用我的。

那麼畫圖就是

原文用柱狀圖展示不同基因在不同處理下的表達量，

但是柱狀圖的信息太少，所以我們用ggplot2做箱線圖。

同樣的方法我們可以做出早花同源異形基因（homeotic gene)，如AP1, APETALA 3(AP3), CAULIFLOWER(CAL)和發育中期和後期基因JAGGED(JAG), WUSCHEL RELATED HOMEOBO 1(WOX1) WOX3的箱線圖：

和文章的結論基本一致，早花同源異形基因在我的流程中只能檢測到CAL, 而發育中期和後期的基因基本上表達量也特別低。不過文章說CAL3這個基因能被穩定被檢測到，實際上表達量與這些基因沒有明顯區別（< 100)，我也不知道他如何得到這個結論的。

文章在比較了不同時期的基因表達情況後，發現

morning-expressed genes of the central oscillator of the circadian clock:LATE ELONGATED HYPOCOTYL (LHY),CRICADOAM CLOCK ASSPCIATED 1 (CCA1)在轉移到長日照下的第一天表達量顯著性提高，而evening-expressed clock genePSEUDO-RESPONSE REGULATOR 5(PRR5)則是顯著性下調。並且evening複合體編碼基因LUX

ARRHYTHMO (LUX),EARLY FLOWERING 4 (ELF4)也和PRR5有相同的趨勢。作者認為這些節律基因表達量的變化是因為日長變化導致的，而不是amplitude改變。