在前面的教程中，我們介紹了使用omicverse完成基本的RNA-seq的分析流程，在本節(jié)教程中侍筛，我們將介紹如何使用omicverse完成加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析WGCNA以及蛋白質(zhì)相互作用PPI分析莽使。
環(huán)境的下載
在這里我們只需要安裝omicverse環(huán)境即可磕瓷，有兩個(gè)方法：
- 一個(gè)是使用conda：
conda install omicverse -c conda-forge - 另一個(gè)是使用pip：
pip install omicverse -i https://pypi.tuna.tsinghua.edu.cn/simple/。-i的意思是指定清華鏡像源沿猜，在國內(nèi)可能會(huì)下載地快一些。
導(dǎo)入環(huán)境
import omicverse as ov
ov.utils.ov_plot_set()
加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)
加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(WGCNA)是一種系統(tǒng)生物學(xué)方法碗脊，用于表征不同樣品之間的基因關(guān)聯(lián)模式啼肩，可用于鑒定高度協(xié)同的基因集，并基于基因集的內(nèi)生性和基因集與表型之間的關(guān)聯(lián)來鑒定候選生物標(biāo)志物基因或治療靶點(diǎn)衙伶。目前引用量已超過15,000祈坠。但Python中完成WGCNA分析相關(guān)的包仍是空白。我們根據(jù)WGCNA的原理矢劲，從底層上復(fù)現(xiàn)了原版WGCNA算法赦拘。
加載數(shù)據(jù)
在這里，我們選擇WGCNA原版的演示數(shù)據(jù)來進(jìn)行分析芬沉，數(shù)據(jù)可以在github上進(jìn)行下載躺同。
import pandas as pd
data=ov.utils.read_csv(filepath_or_buffer='https://raw.githubusercontent.com/Starlitnightly/ov/master/sample/LiverFemale3600.csv',
index_col=0)
data.head()
相關(guān)性矩陣計(jì)算
WGCNA的第一步是計(jì)算基因間的相關(guān)性矩陣，這里我們采用皮爾森系數(shù)的計(jì)算方法丸逸，來完成基因間的直接相關(guān)性矩陣計(jì)算蹋艺。
gene_wgcna=ov.bulk.pyWGCNA(data,save_path='result')
gene_wgcna.calculate_correlation_direct(method='pearson',save=False)
在 pyWGCNA 模塊中，我們需要將直接相關(guān)矩陣轉(zhuǎn)換為間接相關(guān)矩陣來計(jì)算軟閾值黄刚，軟閾值可以幫助我們將原來的相關(guān)網(wǎng)絡(luò)轉(zhuǎn)換為無尺度網(wǎng)絡(luò)
gene_wgcna.calculate_correlation_indirect(save=False)
gene_wgcna.calculate_soft_threshold(save=False)
左邊的垂直坐標(biāo)是無尺度網(wǎng)絡(luò)的評(píng)估指標(biāo)r^2捎谨。R2越接近1，網(wǎng)絡(luò)就越接近無尺度網(wǎng)絡(luò)憔维，通常需要r^2> 0.8或0.9涛救。右側(cè)垂直坐標(biāo)為平均連通度，隨 β 值的增加而減小业扒。將這兩個(gè)圖結(jié)合起來检吆，通常選擇 r^2首次達(dá)到0.8或0.9或更高時(shí)的 β 值。利用 β 值凶赁，我們可以根據(jù)方程將相關(guān)矩陣轉(zhuǎn)換成鄰接矩陣咧栗。
然后我們構(gòu)造拓?fù)渲丿B矩陣。
gene_wgcna.calculate_corr_matrix()
共表達(dá)網(wǎng)絡(luò)分析
在獲得基因間的拓?fù)渲丿B矩陣后虱肄，我們使用動(dòng)態(tài)剪切樹的方式來尋找基因間的模塊致板。在這里，我們使用WGCNA作者發(fā)表的DynamicTree的算法來實(shí)現(xiàn)此功能咏窿。
gene_wgcna.calculate_distance()
gene_wgcna.calculate_geneTree()
gene_wgcna.calculate_dynamicMods()
module=gene_wgcna.calculate_gene_module()
在這里斟或，我們成功地計(jì)算了每個(gè)基因的模塊，共有15個(gè)模塊集嵌，它們的顏色都顯示在圖中萝挤。我們使用.head()查看每個(gè)基因所屬的模塊御毅。
module.head()
我們還可以使用.plot_matrix()來可視化拓?fù)渲丿B矩陣與模塊之間的關(guān)系。
gene_wgcna.plot_matrix()
子模塊分析
有時(shí)候我們對(duì)一個(gè)基因或一個(gè)通路的模塊感興趣怜珍，我們需要提取基因的子模塊進(jìn)行分析和定位端蛆。例如，我們選擇了兩個(gè)模塊6和12作為分析的子模塊
gene_wgcna.get_sub_module([6,12]).shape
我們發(fā)現(xiàn)模塊6和模塊12共有151個(gè)基因酥泛。接下來今豆，我們使用之前構(gòu)建的無尺度網(wǎng)絡(luò)，將閾值設(shè)置為0.95柔袁，為模塊6和模塊12構(gòu)建一個(gè)基因相關(guān)網(wǎng)絡(luò)圖呆躲。
sub_G=gene_wgcna.get_sub_network([6,12],correlation_threshold=0.95)
pyWGCNA 提供了一個(gè)簡(jiǎn)單的可視化函數(shù) plot_sub_network 來可視化我們感興趣的基因相關(guān)性網(wǎng)絡(luò)。
模塊和性狀相關(guān)性分析
除了可以根據(jù)目標(biāo)基因選擇模塊外捶索，我們還可以根據(jù)特定的樣本性狀選擇模塊插掂。我們可以計(jì)算每個(gè)樣本的性狀和模塊之間的相關(guān)性，從而找到與我們感興趣的性狀模塊腥例。
我們首先從 github 中讀取特征矩陣辅甥。特征矩陣形狀必須是以樣本為索引，列為特征院崇。示例名稱必須與前面的原始數(shù)據(jù)的示例名稱一致肆氓。
meta=ov.utils.read_csv(filepath_or_buffer='https://raw.githubusercontent.com/Starlitnightly/ov/master/sample/character.csv',index_col=0)
meta.head()
然后利用分析主成分分析和相關(guān)性，計(jì)算模塊與性狀之間的關(guān)系底瓣。
cor_matrix=gene_wgcna.analysis_meta_correlation(meta)
ax=gene_wgcna.plot_meta_correlation(cor_matrix)
蛋白質(zhì)互作網(wǎng)絡(luò)分析
我們接下來介紹蛋白質(zhì)互作網(wǎng)絡(luò)的分析谢揪，STRING 是一個(gè)已知和預(yù)測(cè)的蛋白質(zhì)-蛋白質(zhì)相互作用的數(shù)據(jù)庫。這種相互作用包括直接的(物理的)和間接的(功能的)關(guān)聯(lián); 它們來源于計(jì)算預(yù)測(cè)捐凭，來源于生物體之間的知識(shí)轉(zhuǎn)移拨扶，以及來自其他(主要的)數(shù)據(jù)庫的相互作用。
在這里茁肠，我們制作了一個(gè)教程患民，使用omicverse來構(gòu)建蛋白質(zhì)-蛋白質(zhì)相互作用網(wǎng)絡(luò)。
演示數(shù)據(jù)
這里我們使用 string-db 的示例數(shù)據(jù)來執(zhí)行分析
酵母中的 FAA4是一種長(zhǎng)鏈脂肪酰輔酶 A 合成酶垦梆，它與其他合成酶以及調(diào)節(jié)劑有關(guān)匹颤。
酵母的NCBI taxonomy Id: 4932
使用omicverse完成蛋白質(zhì)相互作用網(wǎng)絡(luò)分析需要三個(gè)數(shù)據(jù)：蛋白列表
,蛋白類別字典和蛋白顏色字典，顏色字典是繪圖時(shí)的每個(gè)蛋白的顏色托猩，一般與類別字典相同印蓖。在這里我們隨機(jī)生成一個(gè)顏色字典和類別字典。
gene_list=['FAA4','POX1','FAT1','FAS2','FAS1','FAA1','OLE1','YJU3','TGL3','INA1','TGL5']
gene_type_dict=dict(zip(gene_list,['Type1']*5+['Type2']*6))
gene_color_dict=dict(zip(gene_list,['#F7828A']*5+['#9CCCA4']*6))
蛋白質(zhì)相互作用關(guān)系預(yù)測(cè)
omicverse提供了一個(gè)十分簡(jiǎn)單的APIov.bulk.string_interaction京腥，只需要輸入蛋白列表即可完成相互作用關(guān)系的預(yù)測(cè)赦肃。
G_res=ov.bulk.string_interaction(gene_list,4932)
G_res.head()
蛋白質(zhì)相互作用網(wǎng)絡(luò)分析
當(dāng)然，omicverse還有非常漂亮的可視化函數(shù)，來可視化蛋白質(zhì)相互作用網(wǎng)絡(luò)他宛，在這里船侧，我們需要使用pyPPI類來完成上述分析。
#4932代表酵母厅各，人類一般是9606镜撩，小鼠一般是10090
ppi=ov.bulk.pyPPI(gene=gene_list,
gene_type_dict=gene_type_dict,
gene_color_dict=gene_color_dict,
species=4932)
然后我們連接到 string-db 來計(jì)算蛋白質(zhì)-蛋白質(zhì)的相互作用
ppi.interaction_analysis()
我們提供了.plot_network()來可視化蛋白質(zhì)相互作用網(wǎng)絡(luò)，更多的參數(shù)您可以使用help(ov.utils.plot_network)來尋找队塘。
ppi.plot_network()
我們還可以使用ppi.G來獲取蛋白質(zhì)相互作用網(wǎng)絡(luò)琐鲁，該變量的格式為networkx，感興趣的讀者可以自行研究networkx包的相關(guān)分析人灼。
到此，我們的教程就結(jié)束了顾翼，如果你認(rèn)為本教程對(duì)你的研究有幫助投放，不要忘記引用omicverse和WGCNA，最后感謝Jimmy老師對(duì)omicverse的大力支持适贸。
