A. 如何快速從轉錄組數據中篩選目標基因!
如何從海量高通量測序數據中篩選出目標數據?這是困擾大多數老師的一個難題!
我以一個excel的簡單函數為例,演示一下如何從表格中快速篩選感興趣的基因等信息。
函數的名稱是 VLOOKUP函數 ,該函數是Excel表中的一個縱向查找函數,學會該函數之後, 可以方便我們從所有基因的結果中篩選我們關心的基因相關信息,比如:基因的長度、基因在樣品中的表達量、基因的注釋等等信息 。
VLOOKUP函數需要輸入4個值:
1、要查找的值,比如:基因的ID;
2、需要查找的區域;
3、區域中包含返回值的列號,也就是找到相關值之後返回第幾列的信息;
4、精確匹配或者近似匹配,一般我們選擇精確匹配。精確匹配採用0/FALSE、近似匹配採用1/TRUE。
一般來說,我們做完轉錄組測序,都會有一個總表,表裡有所有基因的ID、長度、表達量、差異倍數、注釋信息等等,表格很大,內容很多。
如果我們想提取某些差異基因的基因長度信息,那麼我們該如何操作呢?
我們需要在需要提取長度信息的差異基因表中加上一列gene_length列。
然後插入VLOOKUP函數,按要求輸入4個參數,點擊確定即可。
以上是利用基因ID在總表中查找一列信息,比較簡單。如果我們想查找多列信息該如何操作呢?
方法相似,我們可以在總表中插入deg_gene列,然後去差異基因表中查找基因ID即可,具體操作如下:
到這里,一個簡單的Excel表篩選基因信息的方法就介紹完了,實際上在公司給出轉錄組標准分析之後,很多個性化都可以由自己解決,您需要的僅僅是高手領進門!
更多生物信息課程:
1. 文章越來越難發?是你沒發現新思路,基因家族分析發2-4分文章簡單快速,學習鏈接: 基因家族分析實操課程 、 基因家族文獻思路解讀
2. 轉錄組數據理解不深入?圖表看不懂?點擊鏈接學習深入解讀數據結果文件,學習鏈接: 轉錄組(有參)結果解讀 ; 轉錄組(無參)結果解讀
3. 轉錄組數據深入挖掘技能-WGCNA,提升你的文章檔次,學習鏈接: WGCNA-加權基因共表達網路分析
4. 轉錄組數據怎麼挖掘?學習鏈接: 轉錄組標准分析後的數據挖掘 、 轉錄組文獻解讀
5. 微生物16S/ITS/18S分析原理及結果解讀 、 OTU網路圖繪制 、 cytoscape與網路圖繪制課程
6. 生物信息入門到精通必修基礎課,學習鏈接: linux系統使用 、 perl入門到精通 、 perl語言高級 、 R語言畫圖
7. 醫學相關數據挖掘課程,不用做實驗也能發文章,學習鏈接: TCGA-差異基因分析 、 GEO晶元數據挖掘 、 GSEA富集分析課程 、 TCGA臨床數據生存分析 、 TCGA-轉錄因子分析 、 TCGA-ceRNA調控網路分析
8.其他課程鏈接: 二代測序轉錄組數據自主分析 、 NCBI數據上傳 、 二代測序數據解讀 。
B. 轉錄組WGCNA分析
介紹這個包之前,先要搞清楚這個包能幹啥。(部分內容摘抄自學術咖)
Q1:WGCNA能幹嘛?
A1:能夠將表達模式相似的基因進行聚類,並分析模塊與特定性狀或表型之間的關聯關系。具體一點:1)構建分層聚類樹(hierarchical clustering tree),聚類樹的不同分支代表不同的基因模塊(mole),模塊內基因共表達程度高,而分屬不同模塊的基因共表達程度低。2)探索模塊與特定表型或疾病的關聯關系,最終達到鑒定疾病治療的靶點基因、基因網路的目的。
Q2:WGCNA分析結果中總是提到共表達網路,是什麼?
A2:共表達網路特指利用基因間的表達相關性預測基因間調控關系的方法,WGCNA是共表達網路分析中最有效的方法之一。
Q3:一般說WGCNA的樣品不少於15個,15個樣品考慮重復嗎?
A3:15個樣本這個是包含了生物學重復,比如5個時間點3個重復。
Q4:每個樣本有3個生物學重復,不需要對三個重復的表達量求平均值代表該樣本嗎?
A4:做WGCNA的時候每個樣本是獨立的,三個生物學重復樣本是全部導入做分析,不是取均值再做分析,每個樣本都是獨立的。
Q5:WGCNA裡面一般會提到hubgene,如何確定hubgene?
A5:在WGCNA分析裡面,每個基因都會計算連通性,連通性高的就是hubgene。
那麼根據它能做的事情,再結合具體的數據,那麼我們在做WGCNA之前需要准備的數據有兩個:表達量數據和表型數據。
表達量數據,FPKM矩陣即可。
表型數據,即性狀數據,比如腫瘤的stage、腫瘤的預後等等。可以是質量性狀也可以是數量性狀。
1、安裝包
你可以直接安裝,但是後面會報錯。
看了半天發現,是少了一個impute的包。所以需要重新安裝。
2、導入數據
3、用hclust給所有的樣本建樹。看看不同個體之間的距離,以及有沒有一些具體特別遠的個體。
4、確定最佳的beta值,
畫圖
5(5.1)、構建共表達矩陣(自動構建網路 + 模塊識別)
可視化mole
5(5.2)、構建共表達矩陣(逐漸構建網路 + 模塊識別)
Step_1:Co-expression similarity and adjacency
Step_2:計算拓撲重疊矩陣(TOM)
Step_3:使用TOM(拓撲重疊矩陣)進行聚類,繪制聚類得到的樹形圖。
Step_4:使用dynamic tree cut來識別模塊。
Step_5:將基因表達相似的模塊進行合並
Step_6:保存模塊相關變數,用於後續的分析.需要保存的變數有①模塊的特徵基因②模塊的數字標簽③模塊的顏色標簽④基因的樹形圖。
6、展示模塊之間的相關性
7、可視化基因網路 (TOM plot)
8、模塊和性狀的關聯分析
看完資料之後,性狀關聯分析貌似有兩種處理方法。
第一種:質量性狀。一列subtype但是包含有5種類型的癌症。( https://cloud.tencent.com/developer/article/1516749 )
除了上面的熱圖展現性狀與基因模塊的相關性外。
還可以是條形圖,但是只能是指定某個性狀。
或者自己循環一下批量出圖。
9、感興趣性狀的模塊的相關性分析
C. TCGA 數據分析實戰 —— WGCNA
加權基因共表達網路分析( WGCNA , Weighted gene co-expression network analysis )是一種用來描述不同基因在樣本中的表達關聯模式的系統生物學方法。
通過將表達高度相關的基因聚集成不同的模塊,並探究不同模塊與樣本表型之間的關聯。還可以探究模塊內的關鍵基因的功能,作為潛在的生物標志物或治療靶點進行後續分析
WGCNA 模塊識別演算法大致包含以下幾個步驟:
輸入數據的格式要符合行為樣本,列為基因的矩陣格式,因為計算的是基因之間的相關性,所以數據可以是標准化的表達值或者是 read counts 。
探針集或基因可以通過平均表達量或方差(如中位數或絕對中位差)進行過濾,因為低表達或無變化的基因通常代表噪音。
注意 :並不推薦使用差異基因作為輸入矩陣,通過差異表達基因過濾將會導致一個(或幾個高度相關的)基因聚成一個模塊,同時,也破壞了無標度拓撲的假設,所以通過無標度拓撲擬合來選擇軟閾值的將會失敗。
主要是過濾一些離群或異常的樣本,可以對樣本數據進行聚類,如果存在異常樣本,則其在聚類圖中會顯示出離群現象,可考慮將其剔除。
首先,對基因的表達量進行 0-1 標准化,即
其中, 為樣本方差
然後,使用 pearson 計算基因之間的相關性
兩個基因的共表達相似性表示為
然後將基因之間的相似度轉換為鄰接值,對於非加權網路,計算方式為
其中 為硬閾值,大於等於該閾值表示這兩個基因之間存在連接,而低於閾值則認為兩個基因沒有連接。它們並不能反映共表達信息的連續性質,因此可能導致信息損失。例如,閾值為 0.8 ,那 0.79 是不是應該也有一定的相關性呢?
在介紹軟閾值之前,我們先引出兩個圖論的概念:
度表示為節點所連接的邊的數量
無標度網路具有很好的魯棒性,網路中某些節點的錯誤並不會導致整個網路的癱瘓,具有很多的代償連接。而這一特點,與生物體中的復雜生化網路非常類似,只有少數的基因執行著關鍵性的功能,而大多數的基因執行較為單一的功能。
無標度網路中,節點 d 的度為 k 的概率滿足冪律分布
通過對數變換,變為
從這個公式可以看出,節點的度數與其出現的概率是負相關的,通過計算各個節點的度數 k 與該度數 k 在所有節點度數中的佔比的 pearson 相關性,我們可以得到關於無標度網路的適應系數。該系數越接近 1 則越像無標度網路,越接近 0 則越像隨機網路。
所以,對於加權網路,其鄰接值的計算方式為:
當軟閾值 時,會讓相關系數小的更小,而大的更大。
可以根據適應系數來篩選軟閾值
光有鄰接矩陣是不夠的,基因間的相似性應該要同時體現在其表達和網路拓撲水平,為了能能夠盡可能地最小化噪音和假陽性的影響,因此引入了拓撲重疊矩陣
這個概念的主要表達的是,兩個基因 a 和 b 之間的相關性,不光考慮兩個基因的表達相關性,還需要考慮一些 A 和 B 共有的表達相關基因 u ,如果 u 足夠多,則說明 A 與 B 的網路重疊性強,應該被聚成一類
換個說法,兩個人之間的親密度不僅與他們兩人之間有關,還與他們的共同好友有關,共同好友越多,說明他們兩人之間應該越親密
計算公式為:
其中, 分別為 i 和 j 的度數
表示的是兩個基因的相似性,轉換成距離度量就是 ,並使用該值來進行聚類,並分割模塊
我們以 TCGA 的乳腺癌數據作為示例,來完整的做一遍 WGCNA 分析
先安裝模塊
獲取 50 個樣本的 FPKM 數據, WGCNA 最少需要 15 個樣本, 20 個以上的樣本會更好,樣本越多越好,這里為了方便,我們只挑了 50 個樣本
過濾基因,取絕對中位差 top 5000 的基因
過濾異常樣本
確定軟閾值的時候,需要選擇網路類型,不同的網路類型,其計算鄰接值的方法是不一樣的。
默認為 unsigned
我在 RStudio 中使用 enableWGCNAThreads() 會引發下面的錯誤
所以,我改用了 allowWGCNAThreads() ,就可以運行了
繪制軟閾值曲線
其中橫坐標為軟閾值的梯度,第一幅圖的縱坐標為無標度網路適應系數,越大越好;第二幅圖的縱坐標為節點的平均連通度,越小越好。
查看系統給我們推薦的軟閾值
與我們從圖上看到的結果是一致的,如果出現了異常的值,也就是說在有效的 power 梯度范圍內(無向網路在 power 小於 15 ,有向網路 power 小於 30 ),無法使適應系數的值超過 0.8 ,且平均連接度在 100 以上
可能是由於部分樣品與其他樣品差別較大。這可能是由於批次效應、樣品異質性或實驗條件對表達影響太大等因素造成的。
可以對樣本繪制聚類圖來查看有無異常樣品,如果這確實是由於生物學差異引起的,也可以使用下面的經驗 power 值。
一步法構建網路,我們使用上面推薦的軟閾值 5
查看各模塊的基因數量
可以使用 labels2colors 函數將數值轉換為顏色名稱
使用 plotDendroAndColors 函數來展示各個模塊的層次聚類結果
其中,無法聚類到模塊中的基因會標示為灰色,如果灰色區域較多,可能由於樣本中基因共表達趨勢不明顯,可能需要調整基因過濾的方法。
展示模塊之間的相關性
展示 TOM 矩陣,為了節省時間,我們只使用第一個聚類分支
或者更換一種配色
顏色越深表示基因表達的相關性更高,我們可以看到,模塊內的基因之間具有較高的共表達,而模塊之間的表達相關性較低
將整個網路全部導出成 Cytoscape 輸入文件
保存網路
也可以提取某一模塊的基因
獲取到基因之後,可以進行富集分析找到相關的生物學通路
我們可以分析各網路模塊與樣本表型之間的關系,從而找到與我們感興趣表型相關的模塊。
樣本表型可以是各種指標,比如腫瘤分期分級、已知的分類亞型、葯物響應等,並計算模塊與這些表型之間是否具有顯著相關性
但是模塊是一個矩陣,無法直接計算矩陣和向量之間的相關性,需要轉換為向量之間的相關性。
而 WGCNA 選擇使用 PCA 的方法對數據降維,並將第一主成分定義為 eigengenes ,然後計算 eigengenes 與表型之間的相關性
先獲取並處理臨床數據
計算模塊與 ER 狀態的相關性
如果使用的是其他相關性方法,則可以使用 bicorAndPvalue 函數來計算顯著性
繪制相關性圖
可以看到有些模塊的相關性挺高的,而且也具有顯著性。我們計算出模塊與表型之間相關性之後,可以挑選最相關的那些模塊來進行後續分析。但是,模塊本身可能還包含很多的基因,還需要進一步識別關鍵基因基因。
如何尋找關鍵基因呢?我們可以計算所有基因與模塊之間的相關性,也可以計算基因與表型之間的相關性。如果存在一些基因,既與表型顯著相關又跟某個模塊顯著相關,那麼這些基因可能就是非常重要的關鍵基因了
從上圖中,我們可以看到 paleturquoise 具有較高的相關性,且具有顯著性,我們就來嘗試找找這個模塊的關鍵基因
計算基因與模塊的相關性
再計算基因與表型的相關性
展示模塊內基因與模塊和表型之間的相關性
從圖中我們可以看出,基因與表型的相關性和基因與模塊的相關性還是有一定的線性趨勢的,這說明與表型高度相關的基因,通常也是該表型對應模塊內比較重要的基因。
因此,當我們要選擇關鍵基因時,推薦選取散點圖中右上角部分的基因,即兩個相關性均較大的基因
我們可以導出這個模塊的網路
D. 2021-03-28 WGCNA
加權基因共表達網路分析 (WGCNA, Weighted correlation network analysis)是基於基因的共表達特性進行基因模塊聚類,以探索基因與性狀之間的關聯性,基因模塊與性狀的關聯性,並篩選網路中的核心基因。
相關概念
co-expression network (共表達網路): 一種無方向性的,加權網路,網路的節點代表基因(也可以是蛋白質、代謝產物等),網路的變可以描述基因和基因間共表達程度的高低。為了衡量基因間共表達程度的高低,在計算基因間相關系數(例如皮爾森相關系數)的基礎上,對其進行β次方加權,進而可以強化強相關性節點的關系。
Weighted (加權) :指對相關性值進行冪次運算 。 這種處理方式強化了強相關,弱化了弱相關或負相關,使得相關性數值更符合無標度網路特徵,更具有生物意義。如果沒有合適的 power ,一般是由於部分樣品與其它樣品因為某種原因差別太大導致的,可根據具體問題移除部分樣品或查看後面的經驗值。
Adjacency matrix(鄰接矩陣) :鄰接矩陣有分布在0-1之間的數值組成,是基因和基因之間的加權相關性值構成的矩陣,用來描述節點間相關性強度。
TOM (Topological overlap matrix) :拓撲重疊是通過比較兩個節點和網路中其他節點的加權相關性來定量描述節點間相似性的方法。把鄰接矩陣轉換為拓撲重疊矩陣,以降低噪音和假相關,獲得的新距離矩陣,這個信息可拿來構建網路或繪制TOM圖。
Mole(模塊) :指具有高拓撲重疊相似性的基因集,即高度內連的基因集。共表達模塊是更加非相似性矩陣,利用聚類演算法獲得的。在無向網路中,模塊內是高度 相關 的基因。在有向網路中,模塊內是高度 正相關 的基因。把基因聚類成模塊後,可以對每個模塊進行三個層次的分析:1. 功能富集分析查看其功能特徵是否與研究目的相符;2. 模塊與性狀進行關聯分析,找出與關注性狀相關度最高的模塊;3. 模塊與樣本進行關聯分析,找到樣品特異高表達的模塊。
Mole eigengene (ME) :給定模塊的第一主成分,代表整個模塊的基因表達譜,用來描述模塊在各樣品中的表達模式。
Mole membership (MM) :指給定基因和給定ME之間的相關系數,描述基因屬於一個模塊的可靠性。
Intramolar connectivity (模塊內連通性) :某一個基因的模塊內連通性等同於該基因與模塊內其他基因關聯程度之和,該值越大說明這個基因在模塊中越處於核心位置。
Connectivity (連通性) :類似於網路中 "度"(degree)的概念。每個基因的連連通性是與其相連的基因的邊屬性之和。
Hub gene :關鍵基因 (連接度最多或連接多個模塊的基因)。
Gene significance (GS): 基因顯著性,定義單個基因與外部信息的關聯性,即基因與某個性狀的相關性。
基本分析流程
1.建立關系矩陣:計算兩個基因表達量之間的相關系數,構建成關系矩陣。
2. 建立鄰接矩陣:根據基因表達的相關系數進行加權計算,構建鄰接矩陣。
3. 建立拓撲重疊矩陣:計算節點間的相異程度,將鄰接矩陣轉換為拓撲重疊矩陣。
4.基因模塊識別:基於拓撲鄰接矩陣,進行層級聚類分析,並根據設定標准切分聚類結果,獲得不同的基因模塊,用聚類樹的分枝和不同顏色表示。
5. 核心模塊選擇:根據表型特徵確定核心模塊。
6.核心基因篩選:基於基因連通性篩選核心基因,並圍繞核心基因進行網路構建
WGCNA分析輸入數據
鑒於WGCNA依靠基因的共表達情況進行分析,因此必須要有足夠的樣本數,才能保證相關系數計算的准確性;此外樣本必須包含豐富的變化信息,才能鑒定出有意義的基因模塊。因此WGCNA對於輸入數據有一定的要求:1.不包含生物學重復的獨立樣本組:樣本數>=8;2.包含生物學重復的樣本組:樣本數>=15;3. 輸入數據要求是進行標准化的數據;4. 輸入數據的基因數建議不要超過5000(可以根據變化程度或者表達豐度進行篩選;基因越多,運行時間越長)。
E. WGCNA(轉載)
WGCNA原理及應用
WGCNA介紹:
WGCNA(weighted gene co-expression network analysis,權重基因共表達網路分析)是一種分析多個樣本基因表達模式的分析方法,可將表達模式相似的基因進行聚類,並分析模塊與特定性狀或表型之間的關聯關系,因此在疾病以及其他性狀與基因關聯分析等方面的研究中被廣泛應用。
WGCNA演算法是構建基因共表達網路的常用演算法(詳解: http://www.jianshu.com/p/94b11358b3f3 )。WGCNA演算法首先假定基因網路服從無尺度分布,並定義基因共表達相關矩陣、基因網路形成的鄰接函數,然後計算不同節點的相異系數,並據此構建分層聚類樹(hierarchical clustering tree),該聚類樹的不同分支代表不同的基因模塊(mole),模塊內基因共表達程度高,而分屬不同模塊的基因共表達程度低。最後,探索模塊與特定表型或疾病的關聯關系,最終達到鑒定疾病治療的靶點基因、基因網路的目的。在該方法中mole被定義為一組具有類似表達譜的基因,如果某些基因在一個生理過程或不同組織中總是具有相類似的表達變化,那麼我們有理由認為這些基因在功能上是相關的,可以把他們定義為一個模塊(mole)。這似乎有點類似於進行聚類分析所得到結果,但不同的是,WGCNA的聚類准則具有生物學意義,而非常規的聚類方法(如利用數據間的幾何距離),因此該方法所得出的結果具有更高的可信度。當基因mole被定義出來後,我們可以利用這些結果做很多進一步的工作,如關聯性狀,代謝通路建模,建立基因互作網路等。
WGCNA的用處:
這類處於調控網路中心的基因稱為核心基因(hub gene),這類基因通常是轉錄因子等關鍵的調控因子,是值得我們優先深入分析和挖掘的對象。
在網路中,被調控線連接的基因,其表達模式是相似的。那麼它們潛在有相似的功能。所以,在這個網路中,如果線條一端的基因功能是已知的,那麼就可以預測線條另一端的功能未知的基因也有相似的功能。
下面的問答來自基迪奧,也能加深對WGCNA的理解
問1、調控網路和共表達網路有什麼區別?
答:調控網路是個更廣泛的概念,而共表達網路是調控網路的一種。
理論上我們可以利用各類信息構建調控網路(表達相關性,序列靶向關系、蛋白互作關系),另外調控網路構建的信息既可以來源真實的實驗驗證的關系,也可以來源生物信息的預測。而共表達網路特指利用基因間的表達相關性預測基因間調控關系的方法,而WGCNA又是共表達網路分析中最有效的方法之一。
問2、WGCNA分析適合的生物物種范圍有規定么?
答:沒有限制。對於任何物種中心法則都是存在的,調控關系對於任何物種都是存在的,所以WGCNA沒有物種限定。
問3、同一物種,不同來源的轉錄組數據(比如不同文章/資料來源的),可以放在一起做WGCNA分析嗎?
答:只要樣本間有相似的生物學意義,是可以合並在一起做分析的。但要注意,不同批次之間的樣本是有批次效應的,所以可能會帶來一些誤差,但是是可以放在一起分析的。
問4、相同材料不同處理之間,可以放在一起做WGCNA分析嗎?比如重金屬和鹽鹼處理。
答:可以的。這也正式WGCNA強大的地方,其可以將不同處理的樣本,合並在一起做分析。其他方法則不一定有這么強大的能力,比如做基因表達趨勢分析時,如果樣本涉及到多個處理不同時期的時候,就不好合並分析(或合並後難以解讀)。但WGCNA的方法關注的是調控關系,所以不管是多少個處理組,都可以很好的整合在一起做分析。
問5、不同批次的數據能放一起做WGCNA嗎?
答:可以的。雖然有批次的干擾,但是干擾對WGCNA網路沒有太大影響。因為WGCNA不是做差異分析,而是基因的共表達。因為批次效應理論上不影響相關性。
問6、不同類型的材料,比如親本和F1,適合放一起進行WGCNA么?
答:如果是一個作圖群體,當然親本與F1是可以放在一起分析的,因為你只關心基因的表達模式,所以把親本加進來是沒有問題的。
問7、沒有生物學重復,共3組,每組5個時間點能夠做嗎?
答:理論上有15個樣本,是可以做WGCNA分析的。並且,分析出來的結果對你的研究應該是非常有用的。至少他會比趨勢分析更有意義,更加准確。
問8、一般說WGCNA的樣品不少於15個,15個樣品考慮重復嗎?不同倍性的材料呢?
答:15個樣本這個是包含了生物學重復,比如5個時間點3個重復;在RNA-seq裡面建議不要用不同倍性材料加進來。除非是有參考的多倍體,如果是無參的多倍體,不同倍性之間差異太大,會讓調控網路不準確。所以用單一倍性的材料做調控網路會更加准確。
問9、可以將RNA-seq數據與蛋白組數據,甲基化數據放一起做WGCNA分析?
答:不能與蛋白數據一起分析。因為WGCNA是基於相關系數的演算法。所以最好一起分析的數據變異度是類似的,RNAseq變異非常大,而蛋白的數據變異很小,兩者的變化不在一個數量級上面。所以兩種數據放在一起分析不合理。
但RNA數據可以嘗試跟甲基化數據一起分析。當然我們也建議RNA數據與代謝組數據一起分析,因為代謝組的數據變異也非常大。
問10、表達量和表達的基因數目差異太大的樣品可以一起分析嗎?比如樣品A有2k個gene表達 而樣品B有2w個gene表達了 AB可以一起分析嗎?
答:做WGCNA分析的時候,不能脫離生物學意義,既然要分析調控網路,那麼應該分析有相似生物學意義的一組基因,比如說拿相似組織來一起做分析,比如不應該拿大腦的樣本與腳趾的樣本合並在一起做分析,因為很顯然,這兩個組織沒有關聯。如果兩個樣本之間是有相關聯的生物學意義,哪怕表達的基因數不一樣,或表達模式差異很大,那依然可以放在一起分析;但如果樣本之間完全沒有生物學意義,那麼分析就沒有意義。
問11、實驗設計是case3個時間點(各點都有三個重復),control同樣的3個時間點(每點三個重復),WGCNA怎麼做?3個時間點和case-control兩個因素能同時考慮進來分析嗎?
答:可以的。做WGCNA是更加合理的,因為有兩個梯度的樣本,如果只是做差異分析的話,邏輯可能非常復雜,做WGCNA分析是對樣本特性更好的解析,可以直觀看到基因在六個處理組裡面是怎樣表達的。
問12、可以拿混合樣本分析嗎?比如一個病原細菌跟人類細胞的基因,能說明細菌跟人類細胞基因有調控關系嗎?
答:可以。前提是病原菌有足夠的數據並定量准確,並且這個分析是非常有意義的,最後可以說明這些病原菌可以調控哪些宿主基因。
問13、但是病原宿主混合分析的話,宿主蛋白不能分泌到宿主體內豈不是WGCNA生物學上也沒有意義嗎?
答:依然有意義。即使病原的基因沒有分泌到宿主裡面,但是病原的蛋白是會影響宿主基因的調控的,比如某個細菌感染某個植物,雖然細菌的蛋白不能直接分泌到植物體內,但會影響植物蛋白的分泌。混在一起分析依然是有意義,可以看到植物裡面到底哪個基因對細菌蛋白產生應答作用。
問14、晶元數據兩分類,每組20個樣本,能否每組單獨做WGCNA?
答:可以。WGCNA還有一種重要功能是做兩個網路的比較,比如病人20個樣本做一個調控網路,健康人做一個調控網路,然後兩個網路做比較。
問15、WGCNA可以用來分析lncRNA對下游基因的調控分析嗎?
答:可以。WGCNA網路有利於預測lncRNA的潛在功能。
問16、構建網路是用所有表達基因還是差異基因?
答:這個是具體問題具體分析。如果使用所有的基因分析,會導致運算量非常大。而也不是所有的基因在這個實驗中都有生物學意義,所以我們會提前做一些過濾。
但用於分析的基因不一定是差異表達基因,有時可以用差異表達基因做一個並集,或通過計算變異系數將變異系數低的基因以及低表達的基因去除。但注意,如果你有關心的特定目標基因的話,應該盡量給予保留。
問17、關注某一個pathway上的基因以及調控因子之間的相關性,構建WGCNA網路的時候屬於這個pathway的基因數量太少會不會影響結果呢?
答:這不是問題。在一個調控網路裡面,樣本的某個pathway上,並不是所有基因參與調控(或存在差異性),所以在做WGCNA分析的時候,會做一些過濾,將有變化的基因挑出來再做分析。即分析的是某個pathway上有變化的基因,不需要分析pathway上所有的基因,只需要分析那些變化的基因就夠了。
問18、前期篩選的時候,要選出在所有樣本中變異系數比較大的基因呢?還是直接用差異表達的基因取並集?用基因還是轉錄本,哪個好呢?
答:兩則都可以,我推薦使用變異系數,選擇那些變異較大的基因,來做下面的分析。然後建議用基因不要用轉錄本,因為轉錄本的定量是不準確的。
問19、變異系數一般取多大?
答:具體問題具體分析。例如,沒有特定目標的時候,可以先計算變異系數,將變異系數的百分之前50來做分析,把變異系數偏低的後面一半過濾掉。
問20、輸入數據用FPKM合適嗎?
答:可以。
問21、RNA seq數據是RSEM值怎麼辦?
答:RSEM值原始輸出結果為reads數,如果是RSEM值建議做一個RPKM校正再做分析。
問22、除了RPKM值以外,做WGANA是否還需要其他數據?TCGA數據可否來做WGCNA分析?
答:在做WGCNA分析必須要用表達量數據,但TCGA的數據某些層級沒有表達量數據,沒有表達量數據自然就無法做WGCNA分析。
問23、請問輸入的基因樣本的矩陣的時候,要不要對數據標准化?
答:做WGCNA分析的時候,不需要對數據進行標准化,輸入RPKM值就足以做這個分析。雖然一些文章會做log2處理,但我認為取了LOG2後,會讓一些表達關系沒有那麼豐富。
問24、每個樣本有3個生物學重復,不需要對三個重復的表達量求平均值代表該樣本嗎?
答:注意,做WGCNA的時候每個樣本是獨立的,三個生物學重復樣本是全部導入做分析,不是取均值再做分析,每個樣本都是獨立的。
問25、如果3個生物學重復,做WGCNA的時候是取三個值,還是用cuffdiff處理後取一個值?
答:如果是生物學重復樣本進行調控網路分析,每個樣本獨立使用,而不是取均值。
問26、請問將樣本信息同模塊特徵值進行相關性分析的時候,樣本信息是怎麼處理的呢?比如不同取樣點、不同性別什麼的,這不是數量性狀信息的,這種情況應該怎麼處理呢?
答:樣本的任何信息都可以做模塊相關性分析。比如相關時間點,可以按照先後量化為12134567。又如不同性別,男與女,可以定義為1,-1。任何性狀量化為數字後,都可以進行相關性分析。
問27、怎麼將模塊與性狀對應起來呢有些性狀不好量化,如果直接將模塊與分組對應,如何實現, 不需要量化指標么?
答:首先需要將性狀量化,如果無法將性狀量化,那麼就無法分析。至於分組信息,也可以量化為類似00001111000(1代表一種組別,2代表另一組組別),實現分組信息的數字化。
問28、基因數量為3w左右時,moles數量為多少結果較為理想?怎麼評價聚類效果的好壞?
答:moles數量沒有標准,moles數量無法評估模塊分的好壞,分組是否合理應該看樹的樹形圖,比如樹的分支很清晰就說明模塊式清晰的。moles數量數由生物性狀決定的。比如樣本表達信息很豐富的時候,moles數量會很多;如果樣本的基因表達相對單一,moles數量就會比較少。
問29、我運行例子的時候,得出來基因之間的direction全是undirected,這和前面的幾種關系有什麼區別?
答:WGCNA是一個undirected的方法,它的網路是無方向的,有相關關系但是無方向。
問30、如果做有向網路的構建,您推薦那些方法?
答:很多方法,例如貝葉斯的方法。
問31、非模式物種可以得出基因之間的相互關系類型么?得出的結果也是undirected么?
答:WGCNA是基於表達兩處理的,所以即使是非模式生物,當然也可以他們之間關系,並且關系也是一個無向網路。
問32、選擇幾個表型數據進行結合分析比較好
答:越多越好,看實驗設計。
問33、感染小鼠,5個時間點,3個重復,找不到合適的表型怎麼辦?
答:如果找不到合適表型,可以找某個時間點應答的基因,本身基因的表達趨勢已經有某種生物學意義的。沒有找到合適表型,也可以看變化趨勢。不一定要做表型的相關分析,其他分析也是很有趣的。例如,可以對模塊功能的富集分析,其實都是可以幫助你找到特定模塊的。所以不用糾結於做某個表型的關聯分析。
問34、weight就是tom值嗎?
答:是的。
問35、剪模塊是怎麼做的?是根據TOM劃分嗎?需要自己設定,還是R自動的?
答:剪模塊是R中自動完成的,不需要劃分,但合並的時候你可以設定一個指標,比如差異度是0.25。
問36、看WGCNA說明是用相異矩陣D(D=1-TOM)去做聚類,然後動態剪切?
答:用TOM值來構建矩陣,TOM值就是兩個樣本的相似度,1-TOM值就是兩個樣本的差異度,相似度與差異度可以理解為一個東西,並不矛盾。
問37、模塊特徵值和樣本性狀相關分析的具體方法是?
答:R包用的是計算相關系數的方法。
問38、WGCNA裡面一般會提到hubgene,如何確定hubgene?
答:在WGCNA分析裡面,每個基因都會計算連通性,連通性高的就是hubgene。
問39、在R中安裝「」WGCNA「」說不適合R3.3.1,那適合哪個版本?
答:WGCNA應該是所有版本都適合,如果版本沒有可以考慮降低R軟體的版本,這個對分析沒有影響。因為不同R版本是一樣的。
問40、用STEM分析的時候擬合多少個模型合適?
答:建議不要超過20個。模塊太多不好分析。
參考網站:
http://tiramisutes.github.io/2016/09/14/WGCNA.html
http://www.jianshu.com/p/94b11358b3f3
http://www.omicshare.com/class/home/index/classdetail?id=20
F. wgcna中weight值
你好,首先,wgcna中的weighted(權重)值(也就是題主所說的weigh值)的解釋是:基因之間不僅是相關與否,還保存記錄著它們間的相關性數值,此數值即為基因之間聯系的權重與相關性。然後,wcgna(即weighted gene co-expression network analysis,權重基因共表達網路分析)是適於大樣本分析,分析轉錄組數據的的方法,具體的模型,數據不同,具體的weight值也不同。最後,若想得知所需的weight值,必須要結合具體模型和數據,才可以去計算,分析出其中的weight值,希望可以幫到題主。
G. WGCNA相關重點
WGCNA定義,來源以及發展可閱讀:https://wiki.mbalib.com/wiki/Scale_Free_Network
WGCNA無尺度網路是指在某一復雜的系統中,大部分節點只有少數幾個連結,而某些節點卻擁有與其他節點的大量連結。這些具有大量連結的節點稱為「集散節點」,所擁有的連結可能高達數百、數千甚至數百萬。這一特性說明該網路是無尺度的,因此,凡具有這一特性的網路都是無尺度網路。無尺度網路是指在某一復雜的系統中,大部分節點只有少數幾個連結,而某些節點卻擁有與其他節點的大量連結。這些具有大量連結的節點稱為「集散節點」,所擁有的連結可能高達數百、數千甚至數百萬。這一特性說明該網路是無尺度的,因此,凡具有這一特性的網路都是無尺度網路。
GCNA是加權基因共表達網路分析,旨在分析協同表達的基因模塊,分析基因網路與疾病(表型)之間的關聯,並找出網路中的核心基因。從方法講,WGCNA分為表達量聚類分析和表型關聯。主要包括基因相關系數的計算,基因模塊的確定,共表達網路,模塊與性狀關聯。
R-seq的根本目的是找到差異基因,而在數據歸一化以後所進行的差異基因的尋找有兩種:1 DEseq2差異分析;2 go ontology analysis
第一步 計算基因之間的相關系數 相關系數計算採用相關系數加權值,取相關系數的N次冪,這樣使基因網路呈無尺度網路分布。而閾值的界定決定對相似基因表達的判定。
第二步 構建聚類樹 根據相關系數進行分層聚類,而不同的分枝代表不同的基因模塊,同一基因模塊是表達相似的基因。將幾萬個基因按基因表達相關系數的N次冪分成幾十個模塊。這是信息歸納的過程。
重要含義
鄰接矩陣 就是將頂點的基因和頂點之間的相關系數組成的矩陣稱為鄰接矩陣。一般鄰接矩陣是具體數值,並不是由閾值設定的0/1矩陣。
topilogical matrix
在鄰接矩陣的基礎上再計算一個鄰接矩陣,TOM
WGCNA基本概念
理解WGCNA,需要先理解下面幾個術語和它們在WGCNA中的定義。
共表達網路:定義為加權基因網路。點代表基因,邊代表基因表達相關性。加權是指對相關性值進行冥次運算.(冥次的值也就是軟閾值 (power, pickSoftThreshold這個函數所做的就是確定合適的power))。無向網路的邊屬性計算方式為abs(cor(genex, geney)) ^ power;有向網路的邊屬性計算方式為(1+cor(genex, geney)/2) ^ power; signhybrid的邊屬性計算方式為cor(genex, geney)^power if cor>0 else 0。這種處理方式強化了強相關,弱化了弱相關或負相關,使得相關性數值更符合無標度網路特徵,更具有生物意義。如果沒有合適的power,一般是由於部分樣品與其它樣品因為某種原因差別太大導致的,可根據具體問題移除部分樣品或查看後面的經驗值。
Mole(模塊):高度內連的基因集。在無向網路中,模塊內是高度相關的基因。在有向網路中,模塊內是高度正相關的基因。把基因聚類成模塊後,可以對每個模塊進行三個層次的分析:1. 功能富集分析查看其功能特徵是否與研究目的相符;2. 模塊與性狀進行關聯分析,找出與關注性狀相關度最高的模塊;3. 模塊與樣本進行關聯分析,找到樣品特異高表達的模塊。
基因富集相關文章 去東方,最好用的在線GO富集分析工具;GO、GSEA富集分析一網打進;GSEA富集分析-界面操作。其它關聯後面都會提及。
Connectivity (連接度):類似於網路中 "度"
(degree)的概念。每個基因的連接度是與其相連的基因的邊屬性之和。
Mole eigengene E:
給定模型的第一主成分,代表整個模型的基因表達譜。這個是個很巧妙的梳理,我們之前講過PCA分析的降維作用,之前主要是拿來做可視化,現在用到這個地方,很好的用一個向量代替了一個矩陣,方便後期計算。(降維除了PCA,還可以看看tSNE) (補充:設法將原來變數重新組合成一組新的互相無關的幾個綜合變數,同時根據實際需要從中可以取出幾個較少的綜合變數盡可能多地反映原來變數的信息的統計方法叫做主成分分析或稱主分量分析,也是數學上用來降維的一種方法。)
Intramolar connectivity:
給定基因與給定模型內其他基因的關聯度,判斷基因所屬關系。
Mole membership: 給定基因表達譜與給定模型的eigengene的相關性。
Hub gene: 關鍵基因 (連接度最多或連接多個模塊的基因)。
Adjacency matrix
(鄰接矩陣):基因和基因之間的加權相關性值構成的矩陣。
TOM (Topological overlap
matrix):把鄰接矩陣轉換為拓撲重疊矩陣,以降低噪音和假相關,獲得的新距離矩陣,這個信息可拿來構建網路或繪制TOM圖。
基本分析流程
image
構建基因共表達網路:使用加權的表達相關性。
識別基因集:基於加權相關性,進行層級聚類分析,並根據設定標准切分聚類結果,獲得不同的基因模塊,用聚類樹的分枝和不同顏色表示。
如果有表型信息,計算基因模塊與表型的相關性,鑒定性狀相關的模塊。
研究模型之間的關系,從系統層面查看不同模型的互作網路。
從關鍵模型中選擇感興趣的驅動基因,或根據模型中已知基因的功能推測未知基因的功能。
導出TOM矩陣,繪制相關性圖。
WGCNA包實戰
R包WGCNA是用於計算各種加權關聯分析的功能集合,可用於網路構建,基因篩選,基因簇鑒定,拓撲特徵計算,數據模擬和可視化等。
輸入數據和參數選擇
WGCNA本質是基於相關系數的網路分析方法,適用於多樣品數據模式,一般要求樣本數多於15個。樣本數多於20時效果更好,樣本越多,結果越穩定。
基因表達矩陣:
常規表達矩陣即可,即基因在行,樣品在列,進入分析前做一個轉置。RPKM、FPKM或其它標准化方法影響不大,推薦使用Deseq2的或log2(x+1)對標准化後的數據做個轉換。如果數據來自不同的批次,需要先移除批次效應
(記得上次轉錄組培訓課講過如何操作)。如果數據存在系統偏移,需要做下quantile normalization。
性狀矩陣:用於關聯分析的性狀必須是數值型特徵
(如下面示例中的Height, Weight,
Diameter)。如果是區域或分類變數,需要轉換為0-1矩陣的形式(1表示屬於此組或有此屬性,0表示不屬於此組或無此屬性,如樣品分組信息WT,
KO, OE)。
鏈接:https://www.jianshu.com/p/e9cc3f43441d
下游分析:
得到模塊以後進行 模塊功能富集 計算基因模塊與表型的關系 計算基因與樣本的關系
關鍵挖掘:分析核心基因 利用關系預測基因功能。
實操步驟:
1 數據准備是最復雜的,如果是晶元數據 直接歸一化矩陣即可,而如果是RNAseq數據,那麼用RPKM或者TPM都可以,然後就是樣本的屬性方面信息。
材料准備:需要將正常組織的數據剔除。
2 一般聚類用的是hcust
H. WGCNA分析詳解專題(一)
此次分析詳解專題將講述以下內容,老規矩,如有理解錯誤,還請各位大俠批評指正!
WGCNA適用於什麼分析內容?
表型變數中的分類變數應該如何合理的數值化?
做WGCNA分析我該使用什麼數據,是否需要過濾?
我該選取哪些基因進入分析?是全部的基因還是只用差異表達的基因?
多少樣本量合適呢?怎麼檢測異常(離群樣本)?
如何選取softpower?
如何選取模塊以及模塊中的Hub基因?
WGCNA分析應用(一):發育調控
此次講解應用的文章信息如下:
Title :A novel microglial subset plays a key role in myelinogenesis in developing brain
Published Date :28 September 2017
Published Journal :The EMBO Journal(2017 IF: 10.557)
第一作者 :Agnieszka Wlodarczyk,Department of Neurobiology Research, Institute for Molecular Medicine, University of Southern Denmark(南丹麥大學), Odense, Denmark(丹麥)
1.背景知識
Microglia:小膠質細胞。 中樞神經系統 (central nervous system,CNS) 中的細胞大致分為兩類:神經元(neurons)和神經膠質細胞(glial cells)。小膠質細胞是神經膠質細胞的一種,正常情況下,數量不多,主要分布在大腦、小腦的皮質以及脊髓的灰質中。主要功能:作為中樞神經系統固有的免疫效應細胞,針對刺激,形成活化的小神經膠質細胞,可表達各種抗原,行使抗原遞呈細胞(APC)的功能。
2.數據使用(WGCNA分析使用數據)
數據情況如下:GSE78809(17個樣本)
8個新生兒 Neonates :4個CD11C+和4個CD11C-
6個 EA E(experimental autoimmune encephalomyelitis):3個CD11C+和3個CD11C-
3個成年組 alt :3個CD11C-
解讀 :主要有17個樣本,每一個類別都有大於三個以上的生物學重復,有與大腦發育相關的新生兒組別和成年組別。
3.結果解讀
文章中主要有8個結果,這里我們主要看WGCNA部分的結果,結果3:Distinct gene signatures in microglia subsets ring development and EAE
1.使用的數據:作者使用的是二代測序數據中
所有基因表達的CPM值
WGCNA was applied to the count per million ( CPM ) expression data.
2.圖A: 樣本關系聚類圖 ,這里看到三個組 成年組 , 新生兒組 以及 EAE組 都分開了,並且組內的CD11c+和CD11c-也可以區分開。
疑問點 :有意思的是作者用來做樣本聚類的數值,我在文章找了老久沒有看到圖中橫縱坐標的值是怎麼算的,有知道的可以下方留言討論哈。 一般來說,對樣本進行聚類可以做層次聚類和PCA分析,WGCNA常見的是層次聚類樹。
3.圖B: 模塊聚類樹, 圖的上部分是對基因進行的聚類樹,下面是根據相似性聚成的模塊,文章中總共得到了7個模塊,我們可以在圖E中看到是那幾個模塊以及每個模塊涉及到的功能。
4.圖E: 每個模塊的基因數以及各個模塊的功能 ,灰色模塊是沒有聚類到任何模塊的基因集合。
5.圖C: 表型和模塊相關性圖, 這里可以看出哪些模塊和你關注的表型之間的關系是否顯著
這張圖需要用到一個很重要的表型數據,這里可以看到作者是如何將分類變數數值化的,文章中是這樣描述的:
Six binary variables were generated that were used to calculate the mole trait relationships in which all groups were set to zero with the exception of particular groups of interest:
control (1』s for microglia obtained from healthy control brain)
CD11c (1』s for both EAE CD11c and neonatal CD11c),
EAE (1』s for CD11c negative and microglia obtained from EAE brains),
neonatal (1』s for CD11c negative and microglia obtained from neonatal brains),
CD11c EAE , and CD11c neonatal .
翻譯為表格就是:
對這張圖的解讀很重要,它關繫到了你後面挑選的具體重點分析模塊,以及你如何看這裡面的正相關和負相關,曾經有個小夥伴問我:
這里的負相關算相關么?我看到的大多數文章都是對正相關的結果進行的分析
具體的模塊與表型以及聯合模塊功能的解讀這里就不詳細說了,文章中描述非常詳細,如何將所挖掘到的模塊與發育聯系起來。
4.總結
這篇文章IF在10以上,雖然發表時間比較早了,但是還是值得仔細讀一下的。特別是對結果層面的生物學意義的解讀,很多文章最終結果都只是空泛的說挖掘出了一個biomarker就完了,空洞又無趣
。 作為技術層面的細節,這篇文章里可以看到用於WGCNA分析的目的,樣本數,組內重復樣本數,用來分析的基因,用基因的什麼值,表型數據如何數值化以及對結果如何進行解讀和下游分析。
I. WGCNA分析--提升轉錄組測序文章檔次的利器
現在做轉錄組測序,看看差異基因,做做富集分析,再討論下差異基因功能與自己研究性狀或處理之間的關系,最後加簡單的qPCR驗證,這樣的數據發SCI影響因子越來越低了。必須增加新的分析內容才能有所突破。今天給大家介紹一個能給文章增色的分析內容--基因共表達網路分析(WGCNA),該分析對樣品數有一定要求,建議不少於15個,不過現在測序便宜了,達到這個數量已經不是難事了。下面就給大家介紹兩篇利用WGCNA分析基因共表達網路來提升文章檔次。
文章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
性狀: 大豆籽粒大小
實驗材料
大豆籽粒的大小是一個非常重要的農藝性狀,直接關繫到大豆產量,找到決定大豆籽粒大小的關鍵調控基因對後續的分子育種具有重要意義,因此作者,選取了兩個大豆品種做轉錄組分析,分別是:大籽粒Wandou 28 (V1),小籽粒Peixian Layanghuang (V2),取樣時期為三個時期: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),兩個品種每個樣品三個生物學重復共24個樣品。下圖為種子發育不同時期照片以及籽粒大小差異統計結果:
轉錄組分析結果:
對轉錄組分析結果中每個基因做表達量分析,計算每個基因的表達量FPKM,如果基因的表達量,也就是FPKM值<0.5,認為基因無表達,去除這部分基因。然後,統計每個時期不同品種基因表達量高低的分布圖,大約一半的基因處於低表達水平0.5<=FPKM<=5(下圖A);pca分析發現樣品按照不同發育時期聚類在一起,而不是按照不同品種聚類,說明發育時期是決定基因表達譜的關鍵因素,而性狀的不同引起的轉錄表達差異較小(下圖B),下圖C展示的為不同品種,不同發育時期之間表達基因的韋恩圖,在不同的發育時期都表達的基因還是占絕大多數:
差異基因分析:
差異基因分析,下圖A按相同發育時期,不同的品種之間差異比較,下圖B為不同發育時期之間的差異比較,紅色數字代表上調差異基因數量,黑色代表下調的差異基因數量:
差異基因功能注釋分析,主要針對決定籽粒大小的差異基因的比較,也就是上圖A中的差異基因進行功能分析,挑出一些代表基因,看一下他的功能和表達量,例如,V1S1 vs V2S1差異比較當中,共找到973個差異基因,其中489個基因上調,484個基因下調,上調的代表基因的功能及表達量表格如下圖所示,其中有轉錄因子,植物荷爾蒙(生長素等),脂肪酸代謝,蛋白激酶活性,類黃酮生物合成等功能相關的基因,總之挑選與種子果實等發育生長相關的基因來展示,其他還有好幾個表格,也是關於上圖A中不同時期的上調下調基因的功能注釋表格,展示類似,我這里就不詳細說明了,感興趣的可以查看原文:
不同發育時期差異比較:
不同的發育時期差異基因比較,分別繪制每個發育時期高表達的基因的熱圖,差異基因很多,作者挑選的都是和發育相關,或者和重要農藝性狀相關的差異基因做熱圖,例如轉錄因子相關的基因,荷爾蒙相關的,脂肪酸代謝,澱粉糖代謝等相關的基因。
WGCNA分析找到調控籽粒大小的關鍵hub基因:
首先對所有樣品所有基因的表達量矩陣進行過濾,刪除表達量低的基因(FPKM<0.05),一共有7359個基因用於基因共表達網路構建,總共分析得到12個共表達基因模塊下圖A(聚類樹每一個枝代表一個基因,下面不同的顏色劃分代表基因所處不同的模塊),其中有4個模塊和種子大小相關下圖B,例如,lightyellow模塊,所有的V1的不同時期的樣品與這個模塊高度相關,再例如green模塊,有793個基因,不管是V1樣品,還是V2樣品,這個模塊都與S1相關等等。
4個關鍵模塊基因共表達網路構建發現hub基因:
導出WGCNA共表達網路分析結果,繪制模塊當中基因的表達量熱圖和網路圖,左邊熱圖從上到下分別代表:green mole(A),darkturquoise mole(C),black mole(E),lightyellow mole(G),右邊網路圖分別對應共表達網路,其中紅顏色標記的為連通性較高的hub基因。通過研究這些hub基因的功能發現:這些網路中的關鍵hub基因,包括MYB家族轉錄因子,荷爾蒙(ABA,CK,BA)響應因子,細胞色素P450,BR信號激酶等等,都可能與籽粒的大小相關。
文章2:
題目:
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 天(day after pollination DAP),還測了一下葉片的轉錄組,並取3個生物學重復,共48個樣品。不同發育時期和種子重量差異結果如下:
轉錄組測序結果:
利用轉錄組測序所有基因以及所有樣品的表達矩陣做樣品間的相關性分析和PCA聚類分析,從中可以發現,相同的發育狀態或者組織聚類在一起,說明他們之間具有較強的相關性。
差異基因比較分析:
作者主要比較了相同發育狀態不同品種之間的轉錄組差異比較,差異基因的上下調數量和其中轉錄因子的數量圖a,另外還統計差異基因中不同類型轉錄因子的數量展示圖b,圖c為不同時期差異基因的富集結果,顏色越深說明在該功能上越富集,最後S3時期差異基因在mapman中的Metabolic pathways做了富集分析,可以將差異基因的表達量變化情況展示在通路圖中。
基因共表達網路分析
首先作者將不同的樣品按籽粒大小不同品種分開,分別用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時期)在不同的品種中具有相反的表達,之後作者進一步研究這些模塊裡面基因的相關功能等等:
總結:
上述兩篇文章都是植物當中普通的轉錄組文章,由於添加了WGCNA分析從另一個角度分析與性狀相關的基因,文章的檔次提升不少。想得到WGCNA的分析技能嗎,點擊《 WGCNA視頻教學視頻 》即可觀看:手把手教學包你學會。
更多生物信息課程:
1. 文章越來越難發?是你沒發現新思路,基因家族分析發2-4分文章簡單快速,學習鏈接: 基因家族分析實操課程 、 基因家族文獻思路解讀
2. 轉錄組數據理解不深入?圖表看不懂?點擊鏈接學習深入解讀數據結果文件,學習鏈接: 轉錄組(有參)結果解讀 ; 轉錄組(無參)結果解讀
3. 轉錄組數據深入挖掘技能-WGCNA,提升你的文章檔次,學習鏈接: WGCNA-加權基因共表達網路分析
4. 轉錄組數據怎麼挖掘?學習鏈接: 轉錄組標准分析後的數據挖掘 、 轉錄組文獻解讀
5. 微生物16S/ITS/18S分析原理及結果解讀 、 OTU網路圖繪制 、 cytoscape與網路圖繪制課程
6. 生物信息入門到精通必修基礎課,學習鏈接: linux系統使用 、 perl入門到精通 、 perl語言高級 、 R語言畫圖
7. 醫學相關數據挖掘課程,不用做實驗也能發文章,學習鏈接: TCGA-差異基因分析 、 GEO晶元數據挖掘 、 GSEA富集分析課程 、 TCGA臨床數據生存分析 、 TCGA-轉錄因子分析 、 TCGA-ceRNA調控網路分析
8.其他課程鏈接: 二代測序轉錄組數據自主分析 、 NCBI數據上傳 、 二代測序數據解讀 。
J. 初識WGCNA-基礎知識
WGCNA其譯為 加權基因共表達網路分析 。該分析方法旨在尋找協同表達的基因模塊(mole),並探索基因網路與關注的表型之間的關聯關系,以及網路中的核心基因。
適用於復雜的數據模式, 推薦5組(或者15個樣品)以上的數據 。一般可應用的研究方向有:不同器官或組織類型發育調控、同一組織不同發育調控、非生物脅迫不同時間點應答、病原菌侵染後不同時間點應答。
從方法上來講,WGCNA分為 表達量聚類分析和表型關聯 兩部分,主要包括基因之間相關系數計算、基因模塊的確定、共表達網路、模塊與性狀關聯四個步驟。
第一步計算任意兩個基因之間的相關系數(Person Coefficient)。為了衡量兩個基因是否具有相似表達模式,一般需要設置閾值來篩選,高於閾值的則認為是相似的。但是這樣如果將閾值設為0.8,那麼很難說明0.8和0.79兩個是有顯著差別的。因此, WGCNA分析時採用相關系數加權值,即對基因相關系數取N次冪 ,使得網路中的基因之間的連接服從 無尺度網路分布(scale-freenetworks) ,這種演算法更具生物學意義。
第二步通過基因之間的相關系數構建分層聚類樹, 聚類樹的不同分支代表不同的基因模塊 ,不同顏色代表不同的模塊。基於基因的加權相關系數,將基因按照表達模式進行分類,將模式相似的基因歸為一個模塊。這樣就可以 將幾萬個基因通過基因表達模式被分成了幾十個模塊 ,是一個提取歸納信息的過程。
基因之間不僅僅是相關與否,還記錄著它們的相關性數值,數值就是基因之間的 聯系的權重(相關性)。
模塊(mole):表達模式相似的基因分為一類,這樣的一類基因成為模塊;
Eigengene(eigen + gene):基因和樣本構成的矩陣, https://en.wiktionary.org/wiki/eigengene
鄰近矩陣: 是圖的一種存儲形式,用一個一維數組存放圖中所有頂點數據;用一個二維數組存放頂點間關系(邊或弧)的數據,這個 二維數組 稱為鄰接矩陣;在WGCNA分析裡面指的是基因與基因之間的 相關性系數矩陣。 如果用了閾值來判斷基因相關與否,那麼這個鄰近矩陣就是0/1矩陣,只記錄基因相關與否。但是WGCNA沒有用閾值來卡基因的相關性,而是記錄了所有基因之間的相關性。
WGNA認為基因之間的簡單的相關性不足以計算共表達,所以它利用上面的鄰近矩陣,又計算了一個新的鄰近矩陣。一般來說,TOM就是WGCNA分析的最終結果,後續的只是對TOM的下游注釋。
1.模塊的功能富集
2.模塊與性狀之間的相關性
3.模塊與樣本間的相關系數
1.找到模塊的核心基因
2.利用關系預測基因功能
參考
一文看懂WGCNA 分析(2019更新版)