之前給大家推薦了植物基因家族的分析思路，反響很好，今天還是用文章解讀的方式來剖析轉錄組高階分析思路，希望對大家有所幫助。
隨著轉錄組千元時代的到來，一方面我們可以大膽的去測有價值的樣品而不用擔心成本過高，另一方面審稿專家對相關文章的要求也水漲船高。
想靠著差異基因富集、基因性狀關聯、qPCR驗證這三板斧發好的期刊，難度越來越大了。其實轉錄組數據裡還有很多值得深挖的信息，下面我們就用兩篇發表在《Journal of Experimental Botany》和《The Plant Journal》上的文章來探討下如何通過高階的分析手段來提高文章檔次（文末提供分析方法教學）。
文章1：題目：Identification of regulatory networks and hub genes controlling soybean seed set and size using RNA sequencing analysis
期刊：Journal of Experimental Botany（IF：5.3）
關注性狀：大豆籽粒大小
大豆籽粒的大小是一個非常重要的農藝性狀，直接關係到大豆產量，找到決定大豆籽粒大小的關鍵調控基因對後續的分子育種具有重要意義。
該文作者選取了2個大豆品種做轉錄組測序，分別是：大籽粒Wandou 28 (V1)，小籽粒Peixian Layanghuang (V2)；3個取樣時期分別為：seed set (S1), seed growth (S2), and early seed maturation (S3)；其中前兩個時期的取樣部位分別為：Seed pod with whole seed（S1），Whole seed（S2），S3時期取了兩個部位分別為：Seed coat（S3-1），Seed cotyledon（S3-2），每個樣品3個生物學重複共24個樣品。下圖為種子發育不同時期照片以及籽粒大小差異統計結果：
① 整體表達情況描述：
對每個基因做表達定量，過濾極低表達的基因，統計每個時期不同品種基因表達量高低的分布圖，大約一半的基因處於低表達水平0.5<=FPKM<=5（下圖A）。
PCA分析發現樣品按照不同發育時期聚類在一起，而不是按照不同品種聚類，說明發育時期是決定基因表達譜的關鍵因素，而性狀的不同引起的轉錄表達差異較小（下圖B）。
下圖C展示的為不同品種，不同發育時期之間表達基因的韋恩圖，表明在不同的發育時期都表達的基因還是佔絕大多數：
② 總體差異基因分析：
差異基因分析方面，下圖A按相同發育時期，不同的品種之間比較，下圖B為不同發育時期之間的比較，紅色數字代表上調差異基因數量，黑色代表下調的差異基因數量。
差異基因功能注釋分析，主要針對決定籽粒大小的差異基因，也就是上圖A中的差異基因進行功能分析，挑出一些代表性基因，看一下他的功能和表達量。
例如，V1S1 vs V2S1差異比較當中，共找到973個差異基因（489個上調，484個下調），其中上調基因的功能及表達量表格如下圖所示，其中有轉錄因子，植物荷爾蒙（生長素等），脂肪酸代謝，蛋白激酶活性，類黃酮生物合成等功能相關的基因。總之挑選與種子果實等發育生長相關的基因來展示，其他還有好幾個類似表格，也是關於上圖A中不同時期的上調下調基因的功能注釋表格，這裡就不一一展示了，想詳細了解的可以查看原文。
③ 不同發育時期差異比較：分別繪製每個發育時期高表達基因的熱圖，差異基因很多，作者從中挑選了和發育相關，或者和重要農藝性狀相關的差異基因做熱圖。例如轉錄因子相關的基因，荷爾蒙相關的，脂肪酸代謝，澱粉糖代謝等相關的基因。
WGCNA分析找關鍵hub基因：截止到上面都還是轉錄組常用的分析手段，接下來作者使用了WGCNA（加權基因共表達網絡）分析方法對數據進行了深入挖掘，並成功找到了幾個與籽粒大小相關的關鍵候選基因，成為了本文的亮點，接下來我們看看具體是怎麼做的。
首先對所有樣品所有基因的表達量矩陣進行過濾，刪除表達量低的基因（FPKM<0.05），得到7359個基因用於基因共表達網絡構建。經過分析得到12個共表達基因模塊（見下圖A，聚類樹每一個枝代表一個基因，下面不同的顏色劃分代表基因所處不同的模塊），其中有4個模塊和種子大小相關（見下圖B）。例如lightyellow模塊，所有的V1不同時期的樣品與這個模塊高度相關，再如green模塊，有793個基因，不管是V1樣品，還是V2樣品，這個模塊都與S1相關等等。
導出WGCNA共表達網絡分析結果，繪製4個重要模塊當中基因的表達量熱圖和網絡圖，左邊熱圖從上到下分別代表：green module（A），darkturquoise module（C），black module（E），lightyellow module（G），右邊網絡圖分別對應共表達網絡，其中紅顏色標記的為連通性較高的hub基因。通過研究這些hub基因的功能發現：他們包括MYB家族轉錄因子，荷爾蒙（ABA，CK，BA）響應因子，細胞色素P450，BR信號激酶等等，這些都可能與籽粒的大小相關。
題目：Global transcriptome and co-expression network analyses reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea analysis
期刊：The Plant Journal（IF：5.7）
性狀：鷹嘴豆籽粒大小
這篇文章與上一篇文章思路幾乎一致，只是研究的物種變成了鷹嘴豆。同樣的，也是選取了兩個籽粒大小差異明顯的栽培品種：Himchana 1 (small-seeded) and JGK 3 (large-seeded)，每個品種7個取樣時期S1-S7，分別為授粉後5, 9, 12, 19, 25, 30 and 40 天（DAP, day after pollination），還同時測了葉片的轉錄組，並取3個生物學重複，共48個樣品。不同發育時期和種子重量差異結果如下：
① 整體表達情況描述：
利用轉錄組測序所有基因以及所有樣品的表達矩陣做樣品間的相關性分析和PCA聚類分析，從中可以發現，相同的發育狀態或者組織聚類在一起，說明他們之間具有較強的相關性。
② 差異基因比較分析：
作者主要比較了相同發育狀態不同品種之間的轉錄組差異比較，差異基因的上下調數量和其中轉錄因子的數量（見下圖a），另外還對差異基因中不同類型轉錄因子的數量進行了展示（見下圖b），圖c為不同時期差異基因的富集結果，顏色越深說明在該功能上越富集。最後對S3時期差異基因在mapman中的Metabolic pathways做了富集分析，可以將差異基因的表達量變化情況展示在通路圖中。
WGCNA分析
首先作者將不同的品種按籽粒大小分開，分別進行WGCNA共表達網絡分析，其中在Himchana 1品種中共找到27個模塊(見下圖a)，在JGK 3品種中找到21個模塊(見下圖b)：
進行模塊與樣品之間相關性分析，發現不同發育時期特有的基因模塊（這部分也是分開做）圖中顏色越紅的方框對應的模塊和樣品具有較高的相關性，左邊一半為Himchana 1中模塊與發育時期相關圖，右邊一半為JGK3模塊與發育時期相關結果，然後得到每個樣品中每個時期對應的最相關的模塊，(如下圖)：
結合上一步的分析結果，再來分析兩個品種各自得到的模塊之間的相關性。理論上講，即便是不同品種，其相同發育時期對應的特有模塊應該具有較高的相關性。
例如，通過相關性分析，在JGK 3樣品中左下角黑色模塊與S6發育時期相關，這個模塊與Himchana 1中的darkorange模塊相關，同時darkorange模塊在Himchana 1 中也與S6相關（下圖中紅紫色方框）；同樣的道理其他很多模塊都有這樣的相關性（下圖中紅色方框）。
但是在Himchana 1 中有個orange模塊不與JGK 3中任何一個模塊相關，作者推斷這個特殊的模塊很可能與籽粒大小相關，當然還有其他幾個模塊也有類似的現象。作者進一步研究這些模塊中基因表達情況發現裡面很多基因的表達量（在S3 和 S5時期）在不同的品種中具有相反的表達，之後作者進一步研究這些模塊裡面基因的相關功能等等。
《Journal of Experimental Botany》和《The Plant Journal》都是植物學方面的經典權威期刊，此次解讀的兩篇文章的作者都圍繞研究目的做了多方面的工作。但毫無疑問的是，在文章主體的轉錄組分析部分，WGCNA方法起了重要作用，對數據進行了深入的挖掘，提高了文章的檔次。
如果您也有類似的研究目的，希望對轉錄組數據進行更加深入的利用，那麼WGCNA的方法一定要試一試！
如何實現WGCNA分析經過我們的調研，實現該類型分析有兩種方式：
第一種方式，全套交給公司去做，報價6千到1萬不等，如需調整另外收費；
第二種方式，經過專業的訓練，學習之後自己進行分析工作。已經有部分學術團隊製作了相關的視頻教程，價格在百元級別，能省很多經費，而且自己也多一技在身。
在此給大家推薦一個優秀的相關視頻：組學大講堂製作的《WGCNA-加權基因共表達網絡分析》視頻教程，涵蓋了WGCNA分析的典型內容，並儘量照顧0生信基礎同學，降低了上手難度，一般經過半個月的學習和練習之後就可以掌握WGCNA的分析技能，同時設置了專門的技術交流QQ群在線解答分析中出現的問題。
想用這個方法提升文章檔次的同學趕緊點擊下面連結吧：
http://study.163.com/course/introduction/1005023004.htm?share=1&shareId=1143029466
或者掃描下面的二維碼開始學習之旅！
掃碼開啟學習之旅！！
參考文獻：1. Du J, Wang S, He C, et al. Identification of regulatory networks and hub genes controlling soybean seed set and size using RNA sequencing analysis[J]. Journal of Experimental Botany, 2017, 68(8):1955-1972.
2. Garg R, Singh V K, Rajkumar M S, et al. Global transcriptome and co‐expression network analyses reveal cultivar‐specific molecular signatures associated with seed development and seed size/weight determination in chickpea[J]. Plant Journal, 2017.
延伸閱讀基因家族分析思路解析及學習方法推薦