呦呦切克鬧，WGCNA來一套

劉小澤寫于18.10.10
很早就見過它，但只是不明覺厲，直到昨天才開始認(rèn)真了解

是這樣的：昨天上午看了一篇文獻(xiàn)，講冠心病的基因共表達(dá)分析，簡(jiǎn)而言之就是把得到的幾千個(gè)變化比較大的基因和表型聯(lián)系起來，將他們放入幾十個(gè)模塊中，最后通過比較模塊、表型來確定要研究的那些基因，進(jìn)而可以進(jìn)行功能注釋之類的下游分析。對(duì)于樣本量大的研究（比如芯片、幾十個(gè)樣本的轉(zhuǎn)錄組、GWAS）來說，這么做是最快找到有價(jià)值基因的方法。

文章DOI：10.1186/s12872-016-0217-3
感覺對(duì)自己有幫助，就要拿來學(xué)習(xí)

搜集資料

當(dāng)然，所有pdf資料以及實(shí)戰(zhàn)的教程數(shù)據(jù)都可以在公眾號(hào)回復(fù)“WGCNA”得到

這是一個(gè)全新的領(lǐng)域，但是對(duì)R不陌生，上手會(huì)省時(shí)一些

搜索WGCNA tutorial 第一個(gè)就是啦。最新的教程還是16年更新的，看來作者對(duì)自己的作品比較滿意啊，很少更新了現(xiàn)在

WGCNA官網(wǎng)

踏入新領(lǐng)域，首先就是要了解背景知識(shí)

WGCNA：全稱 Weighted correlation network analysis，中文名是加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

原來做幾個(gè)樣本的轉(zhuǎn)錄組，最后只要得到差異基因，再富集分析一下，就是一個(gè)經(jīng)典的套路；但是樣本數(shù)量如果很多呢？背景基因就會(huì)很多，再用傳統(tǒng)的途徑獲得差異基因，往往會(huì)漏掉許多。WGCNA利用成千上萬個(gè)變化顯著的基因或全部基因的信息來篩選感興趣的基因集，并且聯(lián)系了表型的顯著性關(guān)聯(lián)分析。

這樣做有兩個(gè)好處：一個(gè)是損失的基因少；另一個(gè)是把眾多的基因匯集到幾個(gè)基因集，和表型進(jìn)行關(guān)聯(lián)，就不用做多重假設(shè)檢驗(yàn)

  • Wighted加權(quán)/權(quán)重
    衡量重要程度，也就是不僅要看基因是否相關(guān)，還要看相關(guān)性大小

  • Network 網(wǎng)絡(luò)
    簡(jiǎn)單網(wǎng)絡(luò)：點(diǎn)和線組成，大部分點(diǎn)之間都有線連接
    無尺度網(wǎng)絡(luò)：大部分節(jié)點(diǎn)很少連接，少部分點(diǎn)占據(jù)大量連接?；虻恼{(diào)控關(guān)系就是無尺度網(wǎng)絡(luò)

    簡(jiǎn)單網(wǎng)絡(luò)與無尺度網(wǎng)絡(luò)

  • Correlation相關(guān)
    依據(jù)相關(guān)性得到模塊module ，指的是高度相關(guān)的基因，這些基因在不同組織或者一個(gè)生理過程中有相似的表達(dá)變化。每個(gè)模塊中總會(huì)有佼佼者，也就是“第一成分”，就是說這個(gè)基因和許多基因都有聯(lián)系，可以代表整體模塊。相當(dāng)于北上廣這樣的核心地位

走進(jìn)它的世界

第一步 數(shù)據(jù)輸入與初步處理

一般WGCNA需要15個(gè)樣本以上的數(shù)據(jù)，樣本越多結(jié)果越穩(wěn)定

需要兩樣?xùn)|西：表達(dá)矩陣datExpr、表型矩陣datTraits

  • 表達(dá)矩陣：
    芯片數(shù)據(jù)的歸一化矩陣；
    轉(zhuǎn)錄組的表達(dá)矩陣（一般以RPKM為單位）【一般都是基因在行，樣品在列，后期還需要轉(zhuǎn)置的】

  • 表型矩陣：用于關(guān)聯(lián)分析，必須是數(shù)值型矩陣
    本身就是數(shù)值型：如長度、重量等，直接使用；
    本身不是數(shù)值型，如患病與否這樣的分類變量，需要用0-1表示（0表示沒有該屬性：患病，1表示有該屬性），構(gòu)建0-1矩陣后再分析

    # 舉個(gè)例子
ID    Sick Height Weight 
sam1   1     1   2   
sam2   1     2   7 
sam3   0     10  22 
sam4   0     NA  31  
sam5   0     14  19   

wkdir <- "."
setwd(wkdir)
library(WGCNA)
options(stringsAsFactors = FALSE)
femData = read.csv("LiverFemale3600.csv")
#dim(femData)
names(femData) #看下列名
進(jìn)行轉(zhuǎn)置：設(shè)置初步的表達(dá)矩陣：樣本為行，觀測(cè)為列
datExpr0 = as.data.frame(t(femData[, -c(1:8)]))
names(datExpr0) = femData$substanceBXH
rownames(datExpr0) = names(femData)[-c(1:8)]
檢查基因和樣本是否有太多的缺失值【太多就去掉】
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK # 返回TRUE則繼續(xù)
# （可選）如果存在太多的缺失值
if (!gsg$allOK)
{
  # 把含有缺失值的基因或樣本打印出來
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
  # 去掉那些缺失值
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}
對(duì)樣本進(jìn)行聚類,檢測(cè)異常值
if(T){
  sampleTree = hclust(dist(datExpr0), method = "average")
  pdf(file = "sampleClustering.pdf", width = 12, height = 9)
  par(cex = 0.6)
  par(mar = c(0,4,2,0))
  plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
       cex.axis = 1.5, cex.main = 2)
  # 結(jié)果可以看到，F(xiàn)2_221是一個(gè)異常值
  abline(h = 15, col = "red") #先畫一條輔助線
  dev.off()
}
樣本進(jìn)行聚類
去除異常值，得到過濾后的表達(dá)矩陣
# 把高于15的切除
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust) # 0代表切除的，1代表保留的
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]
開始設(shè)置表型矩陣=》加載表型信息

再次注意：所有的表型信息都要用數(shù)值來表示在各個(gè)樣本中的大小
（比如：體重、長度、肥胖度等都要量化）

traitData = read.csv("ClinicalTraits.csv")
#dim(traitData)
names(traitData)
# 去掉不需要的信息，得到全部樣本的表型矩陣
allTraits = traitData[, -c(31, 16)] #去掉'note'、'comments'
allTraits = allTraits[, c(2, 11:36) ] # 去掉中間的日期之類的
dim(allTraits)
names(allTraits)
得到部分表型矩陣=》表達(dá)矩陣有哪些樣本，表型矩陣就有哪些樣本
femaleSamples = rownames(datExpr)
traitRows = match(femaleSamples, allTraits$Mice)
# 把原來的第一列Mice變成行名
datTraits = allTraits[traitRows, -1]
rownames(datTraits) = allTraits[traitRows, 1]
collectGarbage() #釋放內(nèi)存空間

到這里為止，得到了所有樣本的表達(dá)矩陣和表型矩陣，在進(jìn)一步構(gòu)建網(wǎng)絡(luò)、檢測(cè)模塊之前，先看下表型和樣本的相關(guān)性

針對(duì)去除異常值的表達(dá)矩陣進(jìn)行聚類分析
sampleTree2 = hclust(dist(datExpr), method = "average")
# 用顏色表示表型在各個(gè)樣本的表現(xiàn): 白色表示低，紅色為高，灰色為缺失
traitColors = numbers2colors(datTraits, signed = FALSE)
# 把樣本聚類和表型熱圖繪制在一起
if(T){
  pdf(file = "Sample dendrogram and trait heatmap.pdf", width = 18, height = 10)
  plotDendroAndColors(sampleTree2, traitColors,
                      groupLabels = names(datTraits),
                      main = "Sample dendrogram and trait heatmap")
  dev.off()
}
save(datExpr, datTraits, file = "FemaleLiver-01-dataInput.RData")
樣本聚類和表型熱圖

第二步 網(wǎng)絡(luò)分析

構(gòu)建網(wǎng)絡(luò)、檢測(cè)協(xié)同表達(dá)的基因模塊（非常重要的一步！）

實(shí)現(xiàn)構(gòu)建網(wǎng)絡(luò)有3種方法：(3選1即可)
1.One-step；
2.Step-by-step 支持自定義一些參數(shù)；
3.Block-wise network construction 針對(duì)大型數(shù)據(jù)

加載之前保存的RData

lnames = load(file = "FemaleLiver-01-dataInput.RData")

這里采用一步分析法

篩選軟閾值（soft thresholding power）

原則是使構(gòu)建的網(wǎng)絡(luò)更符合無標(biāo)度網(wǎng)絡(luò)特征

#設(shè)置一系列軟閾值（默認(rèn)1到30）
powers = c(c(1:10), seq(from = 12, to=20, by=2))
#幫助用戶選擇合適的軟閾值，進(jìn)行拓?fù)渚W(wǎng)絡(luò)分析
#需要輸入表達(dá)矩陣、設(shè)置的閾值范圍、運(yùn)行顯示的信息程度(verbose=0不顯示任何信息)
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

sft結(jié)果有兩列，其中powerEstimate返回最佳的軟閾值beta值。這里的結(jié)果是6，表示從6達(dá)到高閾值以后，圖形曲線開始變平，就是說：到達(dá)閾值6時(shí)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)連通的就差不多了，夠本了！具體見下面做出的圖(左一)來理解。實(shí)際使用并不需要知道具體值

畫出結(jié)果 =》橫軸是Soft threshold (power)
(左圖)縱軸是無標(biāo)度網(wǎng)絡(luò)的評(píng)估參數(shù)（相關(guān)系數(shù)的平方),數(shù)值越高，網(wǎng)絡(luò)越符合無標(biāo)度特征 (non-scale)
(右圖)縱軸是基因模塊中所有基因鄰接函數(shù)的均值

if(T){
  pdf(file = "Soft threshold.pdf", width = 18, height = 10)
  par(mfrow = c(1,2))
  cex1 = 0.9
  plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       xlab="Soft Threshold (power)",
       ylab="Scale Free Topology Model Fit,signed R^2",type="n",
       main = paste("Scale independence"))
  #注意這里的-sign(sft$fitIndices[,3])中sign函數(shù)，它把正數(shù)、負(fù)數(shù)分別轉(zhuǎn)為1、-1
  text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       labels=powers,cex=cex1,col="red")
  # 設(shè)置篩選標(biāo)準(zhǔn)h=r^2^=0.9。這里的0.9是個(gè)大概的數(shù)，就是看左圖軟閾值6大概對(duì)應(yīng)的位置
  abline(h=0.90,col="red")
  #看一下Soft threshold與平均連通性
  plot(sft$fitIndices[,1], sft$fitIndices[,5],
       xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
       main = paste("Mean connectivity"))
  text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
  dev.off()
}
得到軟閾值
一步構(gòu)建網(wǎng)絡(luò)[關(guān)鍵一步！]

幾千個(gè)基因組歸類成了幾十個(gè)模塊

計(jì)算基因間的鄰接性，根據(jù)鄰接性計(jì)算基因間的相似性，然后算出基因間的相異性系數(shù)，并因此得到基因間的系統(tǒng)聚類樹
按照混合動(dòng)態(tài)剪切樹的標(biāo)準(zhǔn)，設(shè)置每個(gè)基因模塊最少的基因數(shù)目為30
確定基因模塊后，再次分析，依次計(jì)算每個(gè)模塊的特征向量值
對(duì)模塊進(jìn)行聚類分析，將距離較近的模塊合并為新的模塊
power就是上面計(jì)算得到的軟閾值

net = blockwiseModules(datExpr, power = 6,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "femaleMouseTOM",
                       verbose = 3)
# 顯示模塊數(shù)量以及各自包含的基因數(shù)目
# 0表示未分入任何模塊的基因
# 1是最大的模塊，往后依次降序排列，分別對(duì)應(yīng)各自模塊的基因
table(net$colors)
模塊可視化=》聚類樹

將每個(gè)模塊對(duì)應(yīng)的基因數(shù)轉(zhuǎn)換成顏色單位, 灰色默認(rèn)是無法歸類于任何模塊的那些基因

mergedColors = labels2colors(net$colors)

先畫聚類，后畫顏色

if(T){
  pdf(file = "DendroAndColors.pdf", width = 18, height = 10)
  plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                      "Module colors",
                      dendroLabels = FALSE, hang = 0.03,
                      addGuide = TRUE, guideHang = 0.05)
  dev.off()
}
模塊聚類樹

保存數(shù)據(jù)

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]
save(MEs, moduleLabels, moduleColors, geneTree,
     file = "FemaleLiver-02-networkConstruction-auto.RData")

第三步 模塊聯(lián)系表型信息

看看哪些表型的哪些模塊是自己感興趣的（一般選熱圖顏色深的）

這樣初步表明該模塊中基因可能是有研究價(jià)值的；
下一步，進(jìn)行一個(gè)驗(yàn)證，看看這個(gè)模塊中基因與表型、模塊的相關(guān)性，看看與模塊相關(guān)的基因是不是也與表型相關(guān)

rm(list = ls())
load(file = "FemaleLiver-01-dataInput.RData")
load(file = "FemaleLiver-02-networkConstruction-auto.RData")
模塊關(guān)聯(lián)表型
# 得到基因、樣本數(shù)量
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
# 用color labels重新計(jì)算MEs（Module Eigengenes:模塊的第一主成分）
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p") #（這是重點(diǎn)）計(jì)算ME和表型相關(guān)性
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
對(duì)moduleTraitCor畫熱圖，看結(jié)果挑選自己感興趣的模塊進(jìn)行下游分析
if(T){
  pdf(file = "labeledHeatmap.pdf", width = 18, height = 10)
  # 設(shè)置熱圖上的文字（兩行數(shù)字：第一行是模塊與各種表型的相關(guān)系數(shù)；
  # 第二行是p值）
  # signif 取有效數(shù)字
  textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                     signif(moduleTraitPvalue, 1), ")", sep = "")
  dim(textMatrix) = dim(moduleTraitCor)
  par(mar = c(6, 8.5, 3, 3))
  # 然后對(duì)moduleTraitCor畫熱圖
  labeledHeatmap(Matrix = moduleTraitCor,
                 xLabels = names(datTraits),
                 yLabels = names(MEs),
                 ySymbols = names(MEs),
                 colorLabels = FALSE,
                 colors = greenWhiteRed(50),
                 textMatrix = textMatrix,
                 setStdMargins = FALSE,
                 cex.text = 0.5,
                 zlim = c(-1,1),
                 main = paste("Module-trait relationships"))
  dev.off()
}
模塊關(guān)聯(lián)表型

每種表型都有和它非常相關(guān)的模塊，因此某個(gè)模塊可以作為某個(gè)表型的代表(signature)，對(duì)于非常相關(guān)的模塊，比如weight表型和brown模塊，顏色最深。那么里面的基因是什么？于是進(jìn)行感興趣表型中核心模塊的基因分析

計(jì)算基因與模塊的相關(guān)性矩陣（MM: Module Membership）
# 把各個(gè)module的名字提取出來（從第三個(gè)字符開始），用于一會(huì)重命名
modNames = substring(names(MEs), 3)
# 得到矩陣
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
# 矩陣t檢驗(yàn)
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
# 修改列名
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
計(jì)算基因與表型的相關(guān)性矩陣（GS: Gene Significance）
# 先將感興趣表型weight提取出來，用于計(jì)算矩陣
weight = as.data.frame(datTraits$weight_g)
names(weight) = "weight"
# 得到矩陣
geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"))
# 矩陣t檢驗(yàn)
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
# 修改列名
names(geneTraitSignificance) = paste("GS.", names(weight), sep="")
names(GSPvalue) = paste("p.GS.", names(weight), sep="")
合并兩個(gè)相關(guān)性矩陣,找到和模塊、表型都高度相關(guān)的基因
if(T){
  module = "brown"
  pdf(file = paste0(module,"-MM-GS-scatterplot.pdf"), width = 10, height = 10)
  column = match(module, modNames) #找到目標(biāo)模塊所在列
  moduleGenes = moduleColors==module #找到模塊基因所在行
  par(mfrow = c(1,1))
  verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                     abs(geneTraitSignificance[moduleGenes, 1]),
                     xlab = paste("Module Membership in", module, "module"),
                     ylab = "Gene significance for body weight",
                     main = paste("Module membership vs. gene significance\n"),
                     cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
  dev.off()
}
模塊、表型都高度相關(guān)的基因

圖中看到，MM和GS高度正相關(guān)，因此說明和表型高度相關(guān)的基因，在與表型相關(guān)模塊中也是核心元素。

另外，也可以探索weight表型中其他的模塊，比如第一個(gè)顏色很淺的magenta模塊，它的MM和GS整體負(fù)相關(guān)


負(fù)相關(guān)的基因

第四步 網(wǎng)絡(luò)可視化

【針對(duì)基因】熱圖的方式展示加權(quán)網(wǎng)絡(luò)，每行每列代表一個(gè)基因

還能表示鄰近關(guān)系和拓?fù)渲丿B，淺顏色表示關(guān)系弱，深顏色表示關(guān)系強(qiáng)

基因的聚類分析和模塊顏色畫在頂部和左側(cè)

rm(list = ls())
load(file = "FemaleLiver-01-dataInput.RData")
load(file = "FemaleLiver-02-networkConstruction-auto.RData")
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

選擇400個(gè)基因畫圖

if(T){
  nSelect = 400
  set.seed(10)
  # 計(jì)算拓?fù)渲丿B(TOM: Topological Overlap Matrix)
  # 這個(gè)過程先計(jì)算了鄰接矩陣，后把鄰接矩陣轉(zhuǎn)換為拓?fù)渲丿B矩陣，
  # 降低了噪音和假相關(guān)，獲得距離矩陣dissTOM
  dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6)
  select = sample(nGenes, size = nSelect)
  selectTOM = dissTOM[select, select]
  # 再計(jì)算基因之間的距離樹(對(duì)于基因的子集，需要重新聚類)
  selectTree = hclust(as.dist(selectTOM), method = "average")
  selectColors = moduleColors[select]
  
  pdf(file = paste0("Sub400-netheatmap.pdf"), width = 10, height = 10)
  plotDiss = selectTOM^7
  diag(plotDiss) = NA #將對(duì)角線設(shè)成NA，在圖形中顯示為白色的點(diǎn)，更清晰顯示趨勢(shì)
  TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
  dev.off()
}
基因網(wǎng)絡(luò)可視化

選擇全部基因畫圖(耗時(shí)較久，生成的文件很大)

結(jié)果可以看到各個(gè)區(qū)塊的顏色差異,沿著對(duì)角線的深色區(qū)塊就是模塊Module

if(T){
  pdf(file = "All-netheatmap.pdf", width = 10, height = 10)
  plotTOM = dissTOM^7
  diag(plotTOM) = NA
  TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
  dev.off()
}
【針對(duì)模塊與表型】展示模塊與表型的相關(guān)性

重新計(jì)算模塊的基因樣本相關(guān)矩陣（Eigengenes就是基因和樣本的相關(guān)矩陣）

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
# 將weight表型信息從trait中提取出來
weight = as.data.frame(datTraits$weight_g)
names(weight) = "weight"
# 把weight表型添加到之前計(jì)算的ME矩陣中，并排序
MET = orderMEs(cbind(MEs, weight))

# 畫圖=》meta-modules（模塊的聚類圖加上模塊與表型的熱圖）
# marDendro/marHeatmap 設(shè)置下、左、上、右的邊距
if(T){
  pdf(file = "Eigengene-dengro-heatmap.pdf", width = 10, height = 15)
  par(cex = 0.9)
  plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(4,4,1,2), cex.lab = 0.8, xLabelsAngle
                        = 90)
  dev.off()
}
模塊表型可視化

圖中可以看到，red、blue、brown模塊是高度相關(guān)的，并且它們的相關(guān)性比各自和weight表型的相關(guān)性還大；salmon和weight的相關(guān)性也比較強(qiáng)，但是加入red、blue、brown的meta-modules陣營還不夠格

當(dāng)然分開畫也是可以的

if(T){
  pdf(file = "Eigengene-dendrogram.pdf", width = 10, height = 10)
  par(cex = 1.0)
  plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
                        plotHeatmaps = FALSE)
  dev.off()
  pdf(file = "Eigengene-heatmap.pdf", width = 10, height = 10)
  par(cex = 1.0)
  plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(4,5,2,2),
                        plotDendrograms = FALSE, xLabelsAngle = 90)
  dev.off()
}

第五步 導(dǎo)出網(wǎng)絡(luò)

準(zhǔn)備過程
# 重新計(jì)算拓?fù)渲丿B矩陣
TOM = TOMsimilarityFromExpr(datExpr, power = 6)
# 選擇導(dǎo)出模塊
module = "brown"
# 選擇模塊中基因/探針
probes = names(datExpr)
inModule = (moduleColors==module)
modProbes = probes[inModule]
# 選擇相關(guān)模塊的拓?fù)渲丿B矩陣
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
導(dǎo)出到VisANT
vis = exportNetworkToVisANT(modTOM,
                            file = paste("VisANTInput-", module, ".txt", sep=""),
                            weighted = TRUE,
                            threshold = 0)
感覺基因太多，可以截取部分（比如前30）
nTop = 30
IMConn = softConnectivity(datExpr[, modProbes])
top = (rank(-IMConn) <= nTop)
submodTOM <- modTOM[top, top] #然后用submodTOM代替前面的modTOM導(dǎo)出
導(dǎo)出到Cytoscape
modules = c("brown", "red")
cyt = exportNetworkToCytoscape(modTOM,
                               edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
                               nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
                               weighted = TRUE,
                               threshold = 0.02,
                               nodeNames = modProbes,
                               nodeAttr = moduleColors[inModule])
# 同理，感覺基因太多可以用submodTOM代替modTOM

寫在最后

這只是一個(gè)基礎(chǔ)的流程，為了熟悉熟悉流程。只是針對(duì)一個(gè)大樣本；如果有多個(gè)來源不同的大樣本進(jìn)行分析，還需要參考官網(wǎng)的第二部分教程

得到了模塊中的基因后，就可以進(jìn)行GO、KEGG的富集分析，可以挖掘模塊中關(guān)鍵基因，然后預(yù)測(cè)基因功能

歡迎關(guān)注我們的公眾號(hào)～_～　　
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個(gè)不拽術(shù)語、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容