OneStepWGCNA

雖然題目有點(diǎn)標(biāo)題黨详恼，但這個(gè)插件真的方便，真的感謝這個(gè)插件的作者引几。

工具準(zhǔn)備

TBtools
Rserver: 需要在TBtools的用戶群里下載昧互，一共181M，相當(dāng)于一個(gè)獨(dú)立的R環(huán)境伟桅，具體安裝方法參考Rserver 插件 for TBtools熄浓。用戶群可以在公眾號(hào)：生信藥丸里找到。
OneStepWGCNA插件：可以在TBtools的plugin store下載到勉躺。建議先下載plugin store里的plugin store at high speed缀旁。因?yàn)樵趐lugin store里下載插件太慢了，太慢了盖腕，太慢了赫冬。

數(shù)據(jù)準(zhǔn)備

OneStepWGCNA 的輸入文件是基因表達(dá)矩陣和表型數(shù)據(jù)（或樣品分類數(shù)據(jù)）浓镜。這里用了RED數(shù)據(jù)庫(kù)里的284個(gè)日本晴的不同部位，不同處理的表達(dá)矩陣(FPKM數(shù)據(jù))劲厌。樣品分類數(shù)據(jù)用每個(gè)樣品取的組織膛薛。

表達(dá)矩陣

樣品分類數(shù)據(jù)

第一次跑的時(shí)候遇到的坑

報(bào)錯(cuò)信息

第一次包發(fā)現(xiàn)GenomeInfoDb這個(gè)包沒(méi)下下來(lái)，R文件里也沒(méi)有申明需要這個(gè)包补鼻。我在TBtools的插件目錄下OneStepWGCNA的R腳本里加上了：

if (!require('GenomeInfoDb')) install.packages('GenomeInfoDb')

R文件

再跑一次就成功了
這個(gè)TBtools插件文件夾位置一般在C:\Users\names\ .TBtools.Plugin\OneStepWGCNA

OneStepWGCNA 參數(shù)詳解

readcount or fpkm 是表達(dá)矩陣哄啄，以Tab為分割符，行是樣本风范，列是基因咨跌。
Trait data 是表型文件或分類文件，可以有多個(gè)表型硼婿，可以是連續(xù)性也可以是離散锌半。如果要把多分類性狀和連續(xù)性狀都混到一起，可以把多分類性狀轉(zhuǎn)為one-hot編碼加酵，這個(gè)用下excel就能完成拳喻。

one-hot編碼

Expression data 是表達(dá)矩陣數(shù)據(jù)類型，count數(shù)據(jù)是不能直接用來(lái)做WGCNA的
Normalization method 是表達(dá)矩陣轉(zhuǎn)換的方式猪腕，如果是count數(shù)據(jù)可以選varianceStabilizingTransformation或log10(CPM+1)冗澈，如果是FPKM可以選rawFPKM，不經(jīng)過(guò)轉(zhuǎn)換陋葡，和取對(duì)數(shù)亚亲。建議都試試，看看哪個(gè)效果更好腐缤。
RcCutoff 是用來(lái)除噪聲的捌归，我看原作者的教程里大概意思是count推薦10，F(xiàn)PKM推薦0.2岭粤，當(dāng)然惜索，這個(gè)參數(shù)還是要隨著測(cè)序深度和具體生物學(xué)問(wèn)題而改變的。
samplePerc 和上面的參數(shù)配對(duì)的剃浇，在0-1之間巾兆，如果是0.5，RcCutoff 是0.2 代表想要 50%的樣品的FPKM>0.2虎囚。
RemainGeneNum 是保留多少基因進(jìn)行WGCNA角塑，寫0或小于0就是啥都沒(méi)了，會(huì)報(bào)錯(cuò)淘讥。
mergeCutHeight 是合并模塊的閾值圃伶， 0.25的意思就是把相關(guān)系數(shù)大于0.75的模塊合并。
minModuleSize 是指最少一個(gè)模塊要包含30個(gè)基因。

輸出結(jié)果

Sample_clustering 樣品聚類結(jié)果

樣品聚類結(jié)果

這里如果出現(xiàn)明顯的離群樣本窒朋，記住樣本的名字搀罢，在表達(dá)矩陣和trait矩陣中刪掉對(duì)應(yīng)樣本再跑一次

軟閾值篩選

模塊的鄰接矩陣特征向量相關(guān)性

一共找到了10個(gè)模塊，該結(jié)果主要反應(yīng)模塊間的相關(guān)性拭宁。

模塊劃分

表型與模塊的關(guān)聯(lián)

基因?qū)?yīng)的模塊

最后再嘮叨兩句

我發(fā)現(xiàn)最后結(jié)果里還有很多信息沒(méi)有被存下來(lái)，我想看選幾個(gè)基因畫TOM圖奠货，或者畫個(gè)網(wǎng)絡(luò)圖介褥，或者是圖的字體太小了，看不清想重新畫递惋，該咋辦柔滔？
最粗暴的辦法還是把整個(gè)網(wǎng)絡(luò)保存下來(lái)，想要啥萍虽，自己整睛廊。

修改了下blockwiseModules的參數(shù)，把TOM和net都保存下來(lái)

library(WGCNA)
#把存入的Rdata都讀進(jìn)來(lái)超全，總能用上的 \xk one_traits是我設(shè)的Title name,要改
load("00.one_traits.datatraitbase.Rdata")
load("00.one_traits.net.Rdata")
load("00.one_traits.datatraitbase.Rdata")
load("00one_traits.Modular.Rdata")
load("one_traits.tom-block.1.RData")
Title<-"one_traits"
###TOM圖
moduleColors = labels2colors(net$colors)
TOM <- as.matrix(TOM)
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong 
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function
# 這一部分特別耗時(shí)，行列同時(shí)做層級(jí)聚類
TOMplot(plotTOM, net$dendrograms, moduleColors, 
        main = "Network heatmap plot, all genes")

###導(dǎo)出Cytoscape 圖格式
probes = colnames(datExpr)
dimnames(TOM) <- list(probes, probes)
# Export the network into edge and node list files Cytoscape can read
# threshold 默認(rèn)為0.5, 可以根據(jù)自己的需要調(diào)整王财，也可以都導(dǎo)出后在
# cytoscape中再調(diào)整
cyt = exportNetworkToCytoscape(TOM,
             edgeFile = paste(Title, ".edges.txt", sep=""),
             nodeFile = paste(Title, ".nodes.txt", sep=""),
             weighted = TRUE, threshold = 0.5,
             nodeNames = probes, nodeAttr = moduleColors)

## 模塊內(nèi)基因與表型數(shù)據(jù)關(guān)聯(lián)
#關(guān)聯(lián)樣品性狀的二元變量時(shí)卵迂，設(shè)置
robustY = ifelse(corType=="pearson",T,F)
if (corType=="pearsoon") {
  geneModuleMembership = as.data.frame(cor(datExpr, MEs_col, use = "p"))
  MMPvalue = as.data.frame(corPvalueStudent(
             as.matrix(geneModuleMembership), nSamples))
} else {
  geneModuleMembershipA = bicorAndPvalue(datExpr, MEs_col, robustY=robustY)
  geneModuleMembership = geneModuleMembershipA$bicor
  MMPvalue   = geneModuleMembershipA$p
}
# 計(jì)算性狀與基因的相關(guān)性矩陣
## 只有連續(xù)型性狀才能進(jìn)行計(jì)算，如果是離散變量绒净，在構(gòu)建樣品表時(shí)就轉(zhuǎn)為0-1矩陣见咒。
traitData<-read.table(traitsfile,header=T,sep='\t') #traitsfile 表型文件
if (corType=="pearsoon") {
  geneTraitCor = as.data.frame(cor(datExpr, traitData, use = "p"))
  geneTraitP = as.data.frame(corPvalueStudent(
             as.matrix(geneTraitCor), nSamples))
} else {
  geneTraitCorA = bicorAndPvalue(datExpr, traitData, robustY=robustY)
  geneTraitCor = as.data.frame(geneTraitCorA$bicor)
  geneTraitP   = as.data.frame(geneTraitCorA$p)
}

## Warning in bicor(x, y, use = use, ...): bicor: zero MAD in variable 'y'.
## Pearson correlation was used for individual columns with zero (or missing)
## MAD.

# 最后把兩個(gè)相關(guān)性矩陣聯(lián)合起來(lái),指定感興趣模塊進(jìn)行分析
module = "brown"#感興趣的模塊
pheno = "Root" #感興趣的性狀
modNames = substring(colnames(MEs_col), 3)
# 獲取關(guān)注的列
module_column = match(module, modNames)
pheno_column = match(pheno,colnames(traitData))
# 獲取模塊內(nèi)的基因
moduleGenes = moduleColors == module

# 與性狀高度相關(guān)的基因，也是與性狀相關(guān)的模型的關(guān)鍵基因
verboseScatterplot(abs(geneModuleMembership[moduleGenes, module_column]),
                   abs(geneTraitCor[moduleGenes, pheno_column]),
                   xlab = paste("Module Membership in", module, "module"),
                   ylab = paste("Gene significance for", pheno),
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

10000個(gè)基因的TOM

棕色模塊與根的相關(guān)性