[#R] WGCNA 保守性分析

最新推荐文章于 2024-08-25 20:50:23 发布
最新推荐文章于 2024-08-25 20:50:23 发布
阅读量305 收藏 3
点赞数 7
文章标签： r语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/parkzx/article/details/140316661
版权

结果示例：

参考：

https://www.jianshu.com/p/9a03a35cabf3

https://www.jianshu.com/p/e9e2ef9ea693

输入文件示例：

#人和小鼠 ID 转化

Combining species datasets

#hsa.exp.2.xls#

#mmu17.exp.2.xls#

#mmu80.exp.2.xls##

代码：

library(WGCNA)
setwd("E:/GEO")

hsa = read.delim('hsa.exp.2.xls', row.names = 1, sep = '\t', check.names = FALSE)
mmu17=read.delim('mmu17.exp.2.xls', row.names = 1, sep = '\t', check.names = FALSE)
mmu80=read.delim('mmu80.exp.2.xls', row.names = 1, sep = '\t', check.names = FALSE)

hsa2 =hsa[ , -which(colnames(hsa) %in% c("Cntrl1","Cntrl2"))]
#sampleTree=hclust(dist(t(hsa)),method = "average")
#par(cex=0.5)
#plot(sampleTree)

#hsa2=hsa[ , -which(colnames(hsa) %in% c("C5","",""))]  #C5离群，去除

#sampleTree=hclust(dist(t(mmu17)),method = "average")
#par(cex=0.5)
#plot(sampleTree)

mmu17_2=mmu17 #没有


mmu80_2=mmu80[ , -which(colnames(mmu80) %in% c("GSM1904766","GSM1904768","GSM1904769","GSM1904771","GSM1904772","GSM1904773",
"GSM1904775","GSM1904777","GSM1904778","GSM1904780","GSM1904781","GSM1904782","GSM1904784","GSM1904785","GSM1904786","GSM1904788","GSM1904790",
"GSM1904791","GSM1904792","GSM1904793","GSM1904800","GSM1904801","GSM1904803","GSM1904804","GSM1904806","GSM1904814",
"GSM1904815","GSM1904816","GSM1904817","GSM1904818","GSM1904822","GSM1904823","GSM1904825","GSM1904826","GSM1904827","GSM1904829",
"GSM1904833","GSM1904834","GSM1904838","GSM1904839"))]

powers = c(c(1:10), seq(from = 12, to=50, by=2))
sft = pickSoftThreshold(t(mmu80_2), powerVector = powers, verbose = 5)
cex1 = 0.5
##画出
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")) +
  text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       labels=powers,cex=cex1,col="red")+
  abline(h=0.90,col="red")
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")
hsa2_power=18
mmu17_2_power=26
mmu80_2_power=10

hsa_net = blockwiseModules(
  ##这个矩阵行是样本，列是基因
  t(hsa2), 
  power = 18,
  TOMType = "unsigned", minModuleSize = 40,
  reassignThreshold = 0, mergeCutHeight = 0.25,
  numericLabels = TRUE, pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  #saveTOMFileBase = "MyTOM",
  verbose = 3)
hsa_mergedColors = labels2colors(hsa_net$colors)
pdf("hsa_module.pdf",width = 10, height = 5)
plotDendroAndColors(hsa_net$dendrograms[[1]], hsa_mergedColors[hsa_net$blockGenes[[1]]], "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()


mmu17_net = blockwiseModules(
  ##这个矩阵行是样本，列是基因
  t(mmu17_2), 
  power = 26,
  TOMType = "unsigned", minModuleSize = 40,
  reassignThreshold = 0, mergeCutHeight = 0.25,
  numericLabels = TRUE, pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  #saveTOMFileBase = "MyTOM",
  verbose = 3)
mmu17_mergedColors =labels2colors(mmu17_net$colors)
pdf("mmu17_module.pdf",width = 10, height = 5)
plotDendroAndColors(mmu17_net$dendrograms[[1]], mmu17_mergedColors[mmu17_net$blockGenes[[1]]], "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

mmu80_net = blockwiseModules(
  ##这个矩阵行是样本，列是基因
  t(mmu80_2), 
  power = 10,
  TOMType = "unsigned", minModuleSize = 40,
  reassignThreshold = 0, mergeCutHeight = 0.25,
  numericLabels = TRUE, pamRespectsDendro = FALSE,
  saveTOMs = TRUE,
  #saveTOMFileBase = "MyTOM",
  verbose = 3)
mmu80_mergedColors = labels2colors(mmu80_net$colors)
pdf("mmu80_module.pdf",width = 10, height = 5)
plotDendroAndColors(mmu80_net$dendrograms[[1]], mmu80_mergedColors[mmu80_net$blockGenes[[1]]], "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

#table(hsa_net$colors)
#labels2colors(net$colors)
hsa_mergedColors = labels2colors(hsa_net$colors)
mmu17_mergedColors =labels2colors(mmu17_net$colors)
mmu80_mergedColors = labels2colors(mmu80_net$colors)



#hsa+mmu17
nSets = 2
multiExpr = list()
hsa2 = subset(hsa2,rowSums(hsa2) >1) 
#write.table(hsa2,file="hsa.tmp,xls",sep = "\t")
mmu17_2=subset(mmu17_2,rowSums(mmu17_2) >1)
multiExpr[[1]] = list(data = t(hsa2))
multiExpr[[2]] = list(data = t(mmu17_2))
setLabels = c("Human", "Mouse17")
names(multiExpr) = setLabels
lapply(multiExpr, lapply, dim)

multiColor = list(Human = hsa_mergedColors,Mouse17=mmu17_mergedColors)
#colorList = list(hsa_mergedColors, mmu17_mergedColors)
#names(colorList) = setLabels

mp = modulePreservation(multiExpr, multiColor,
                        referenceNetworks = c(1:2),
                        loadPermutedStatistics = FALSE,
                        nPermutations = 200,
                        verbose = 3)
ref = 1
test = 2
statsObs = cbind(mp$quality$observed[[ref]][[test]][, -1], mp$preservation$observed[[ref]][[test]][, -1])
statsZ = cbind(mp$quality$Z[[ref]][[test]][, -1], mp$preservation$Z[[ref]][[test]][, -1])
print(cbind(statsObs[, c("medianRank.pres", "medianRank.qual")],
            signif(statsZ[, c("Zsummary.pres", "Zsummary.qual")], 2)) )
modColors = rownames(mp$preservation$observed[[ref]][[test]])
moduleSizes = mp$preservation$Z[[ref]][[test]][, 1]
row.names(statsZ[statsZ$Zsummary.pres<10,])

#去掉Z<10的module
#%in%不在这里的基因
plotMods = !(modColors %in% row.names(statsZ[statsZ$Zsummary.pres<-10,]))
#去掉了Z<10的基因
text = modColors[plotMods]
plotData = cbind(mp$preservation$observed[[ref]][[test]][, 2], mp$preservation$Z[[ref]][[test]][, 2])
##preservation可视化
mains = c("Preservation Median rank", "Preservation Zsummary")
##新开一个画图窗口
sizeGrWindow(10, 5)
pdf("hsa_mmu17_preservation.pdf",width = 20, height = 10)
##一行两列
par(mfrow = c(1,2))
##到四边的距离
par(mar = c(4.5,4.5,2.5,1))
for (p in 1:2){
  min = min(plotData[, p], na.rm = TRUE);
  max = max(plotData[, p], na.rm = TRUE);
  # Adjust ploting ranges appropriately
  if (p==2){
    if (min > -max/10) min = -max/10
    ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
  } else
    ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
  #bg 颜色 pch 圆圈的种类
  plot(moduleSizes[plotMods], plotData[plotMods, p], col = 1, bg = modColors[plotMods], pch = 21,
       main = mains[p],
       ##圆圈的大小
       cex = 2.4,
       ylab = mains[p], xlab = "Module size", log = "x",
       ylim = ylim,
       xlim = c(10, 2000), cex.lab = 1.2, cex.axis = 1.2, cex.main =1.4)
  ##贴上标签
  labelPoints(moduleSizes[plotMods], plotData[plotMods, p], text, cex = 1, offs = 0.08);
  # For Zsummary, add threshold lines
  if (p==2){
    abline(h=0)
    abline(h=2, col = "blue", lty = 2)
    abline(h=10, col = "darkgreen", lty = 2)
  }
}
dev.off()
print(signif(statsZ[, "Zsummary.pres", drop = FALSE],2))
overlap = cbind(mp$accuracy$observedCounts[[1]][[2]], mp$accuracy$observedFisherPvalues[[1]][[2]])
print(overlap)

#has+mmu80
nSets = 2
multiExpr = list()
hsa2 <- subset(hsa2,rowSums(hsa2) >1)

mmu80_2=subset(mmu80_2,rowSums(mmu80_2) >1)
multiExpr[[1]] = list(data = t(hsa2))
multiExpr[[2]] = list(data = t(mmu80_2))
setLabels = c("Human", "Mouse80")
names(multiExpr) = setLabels
lapply(multiExpr, lapply, dim)

multiColor = list(Human = hsa_mergedColors,Mouse80=mmu80_mergedColors)
#colorList = list(hsa_mergedColors, mmu17_mergedColors)
#names(colorList) = setLabels

mp_hsa_mm80 = modulePreservation(multiExpr, multiColor,
                        referenceNetworks = c(1:2),
                        loadPermutedStatistics = FALSE,
                        nPermutations = 200,
                        verbose = 3)

ref = 1
test = 2
statsObs = cbind(mp_hsa_mm80$quality$observed[[ref]][[test]][, -1], mp_hsa_mm80$preservation$observed[[ref]][[test]][, -1])
statsZ = cbind(mp_hsa_mm80$quality$Z[[ref]][[test]][, -1], mp_hsa_mm80$preservation$Z[[ref]][[test]][, -1])
print(cbind(statsObs[, c("medianRank.pres", "medianRank.qual")],
            signif(statsZ[, c("Zsummary.pres", "Zsummary.qual")], 2)) )
modColors = rownames(mp_hsa_mm80$preservation$observed[[ref]][[test]])
moduleSizes = mp_hsa_mm80$preservation$Z[[ref]][[test]][, 1]
#row.names(statsZ[statsZ$Zsummary.pres<10,])

#去掉Z<10的module
#%in%不在这里的基因
plotMods = !(modColors %in% row.names(statsZ[statsZ$Zsummary.pres<-10,]))
#去掉了Z<10的基因
text2 = modColors[plotMods]
plotData = cbind(mp_hsa_mm80$preservation$observed[[ref]][[test]][, 2], mp_hsa_mm80$preservation$Z[[ref]][[test]][, 2])
##preservation可视化
mains = c("Preservation Median rank", "Preservation Zsummary")
##新开一个画图窗口
sizeGrWindow(10, 5)
pdf("hsa_mmu80_preservation.pdf",width = 20, height = 10)
##一行两列
par(mfrow = c(1,2))
##到四边的距离
par(mar = c(4.5,4.5,2.5,1))
for (p in 1:2){
  min = min(plotData[, p], na.rm = TRUE);
  max = max(plotData[, p], na.rm = TRUE);
  # Adjust ploting ranges appropriately
  if (p==2){
    if (min > -max/10) min = -max/10
    ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
  } else
    ylim = c(min - 0.1 * (max-min), max + 0.1 * (max-min))
  #bg 颜色 pch 圆圈的种类
  plot(moduleSizes[plotMods], plotData[plotMods, p], col = 1, bg = modColors[plotMods], pch = 21,
       main = mains[p],
       ##圆圈的大小
       cex = 2.4,
       ylab = mains[p], xlab = "Module size", log = "x",
       ylim = ylim,
       xlim = c(10, 2000), cex.lab = 1.2, cex.axis = 1.2, cex.main =1.4)
  ##贴上标签
  labelPoints(moduleSizes[plotMods], plotData[plotMods, p], text2, cex = 1, offs = 0.08);
  # For Zsummary, add threshold lines
  if (p==2){
    abline(h=0)
    abline(h=2, col = "blue", lty = 2)
    abline(h=10, col = "darkgreen", lty = 2)
  }
}
dev.off()
print(signif(statsZ[, "Zsummary.pres", drop = FALSE],2))
overlap = cbind(mp_hsa_mm80$accuracy$observedCounts[[1]][[2]], mp_hsa_mm80$accuracy$observedFisherPvalues[[1]][[2]])
print(overlap)

table(hsa_net$colors)
table(mmu17_net$colors)
table(mmu80_net$colors)

for(i in 1:length(text)){
  y=hsa2[which(hsa_mergedColors==text[i]),]
  write.csv(y,paste(paste("hsa_mm17/hsa.",text[i],sep = ""),"csv",sep = "."),quote=F)
}

for(i in 1:length(text)){
  y=mmu17_2[which(mmu17_mergedColors==text[i]),]
  write.csv(y,paste(paste("hsa_mm17/mmu17.",text[i],sep = ""),"csv",sep = "."),quote=F)
}


for(i in 1:length(text2)){
  y=hsa2[which(hsa_mergedColors==text2[i]),]
  write.csv(y,paste(paste("hsa_mm80/hsa.",text2[i],sep = ""),"csv",sep = "."),quote=F)
}

for(i in 1:length(text2)){
  y=mmu80_2[which(mmu80_mergedColors==text2[i]),]
  write.csv(y,paste(paste("hsa_mm80/mmu80.",text2[i],sep = ""),"csv",sep = "."),quote=F)
}

我的飞鸿踏雪泥
关注
  • 7
    点赞
  • 踩
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
氨基酸保守性.zip:利用保守指数分析氨基酸保守性-matlab开发
05-30
考虑到氨基酸之间的生化相似性，我们分析了多序列比对 (MSA) 中每个位置的保守程度。 在我们的实验中，我们将守恒指数 (CI) 与其他最先进的守恒方法（即香农熵、Jensen-Shannon 多样性得分 [2]、相对熵 [2]、冯...
AnoConEly:昆虫Dscam基因内的保守性分析
05-11
昆虫Dscam基因内的保守性分析Watly（at）knox.edu的Liam Ely Contact的DSCAM研究 从NCBI检索到的序列数据，可通过NCBI的BLASTN / BLAST +找到。 激发这项研究的论文： : 非常感谢Matthew Jones-Rhoades博士和Ole ...
用R语言做WGCNA分析全步骤二（附代码解读）【转载】
qq_43138237的博客
04-11 5771
代码逐句分析一、文章来源二、 关联模块与临床特征1.量化module-trait(模块-特征)关系二级目录三级目录 一、文章来源 初学WGCNA，觉得博主的代码写的很不错，但其中很多代码在第一遍看的时候有很多地方不理解，后查阅了很多资料，终于看明白了，于是写了一篇笔记，记录自己的学习心得，有不准确的地方，还望各位大佬们不吝赐教~ 原文链接：WGCNA简明指南 二、 关联模块与临床特征 1.量化module-trait(模块-特征)关系 在这个分析中，我们希望识别与临床特征显著相关的模块。由于我们已经有了每个
线下开始 | 高级转录组分析和R数据可视化
悟道西方
03-02 463
福利公告：为了响应学员的学习需求，经过易生信培训团队的讨论筹备，现安排《高级转录组分析和R数据可视化》于2023年3月10-12线上/线下课程 (线上课是通过腾讯会议实时直播线下课，实时互动，并录制有视频回放，无限期观看)。报名参加线上直播课的老师可在365天内选择参加同课程的一次线下课。期待和大家的线上线下相识。转录组分析是目前应用最广的高通量测序分析技术之一。常见设计是不同样品之间比较，寻找...
转录组分析的正确姿势你了解了吗？
刘永鑫的博客——宏基因组公众号
08-27 950
转录组分析是目前应用最广的高通量测序分析技术之一。常见设计是不同样品之间比较，寻找差异基因、标志基因、协同变化基因、差异剪接和新转录本，并进行结果可视化、功能注释和网络分析等。转录组的测序...
经典入门 | 高级转录组分析和R数据可视化 （2024.9）
刘永鑫的博客——宏基因组公众号
06-03 396
福利公告：为了响应学员的学习需求，经过易生信培训团队的讨论筹备，现安排《高级转录组分析和R数据可视化》于2024年9月20-22线上/线下课程 (线上课是通过腾讯会议实时直播线下课，实时互动，并录制有视频回放，无限期观看)。报名参加线上直播课的老师可在365天内选择参加同课程的一次线下课。期待和大家的线上线下相识。相关课程转录组线上/线下开课时间：2024/9/20-22扩增子线上/线下开课时间...
高级转录组分析和R数据可视化技术研讨会(2023.2)
刘永鑫的博客——宏基因组公众号
11-16 198
福利公告：为了响应学员的学习需求，经过易生信培训团队的讨论筹备，现安排《高级转录组分析和R数据可视化》于2023年2月17线上/线下课程 (线上课是通过腾讯会议实时直播线下课，实时互动，并录制有视频回放，无限期观看)。报名参加线上直播课的老师可在365天内选择参加同课程的一次线下课。期待和大家的线上线下相识。扩增子线上/线下开课时间：2023/2/3宏基因组线上/线下开课时间：2022/11/1...
相约2023 | 第 20 期高级转录组分析和R数据可视化火热报名中！！！
悟道西方
11-24 176
福利公告：为了响应学员的学习需求，经过易生信培训团队的讨论筹备，现安排《高级转录组分析和R数据可视化》于2023年02月17-19 线上/线下课程 (线上课是通过腾讯会议实时直播线下课，实时互动，并录制有视频回放，无限期观看)。报名参加线上直播课的老师可在365天内选择参加同课程的一次线下课。期待和大家的线上线下相识。临床基因组学线上/线下开课时间：2022/11/4-2022/11/6宏基因组...
生信分析—差异分析,聚类分析，相关性分析
anyinglengtong的博客
11-18 5435
生信中提到的二代数据通常指的是第二代测序数据，即在基因组学研究中使用的第二代高通量测序技术生成的数据。第二代测序技术包括Illumina/Solexa、454/Roche、Ion Torrent等平台，它们通常以高通量、低成本、高准确性和较短的读长为特点。这些技术的出现使得大规模基因组测序成为可能，对于遗传学、演化生物学、临床医学等领域的研究起到了重要的推动作用。表型：生物个体可观测的性状。基因型：是指某一生物个体全部基因组合的总称。它反映生物体的遗传构成，即从双亲获得的全部基因的总和。
生信入门转录组和可视化学习捷径
悟道西方
08-29 1048
转录组分析是目前应用最广的高通量测序分析技术之一。常见设计是不同样品之间比较，寻找差异基因、标志基因、协同变化基因、差异剪接和新转录本，并进行结果可视化、功能注释和网络分析等。转录组的测序...
微生物组-宏基因组分析专题研讨会(2022.8)
刘永鑫的博客——宏基因组公众号
06-02 767
福利公告：为了响应学员的学习需求，经过易生信培训团队的讨论筹备，现决定安排扩增子16S分析、宏基因组、Python课程和转录组的线上直播课。报名参加线上直播课的老师可在1年内选择参加同课程的一次线下课。期待和大家的线上线下会晤。目前可以通报的信息：临床基因组学线上/线下开课时间：2022/7/8-2022/7/10扩增子线上/线下开课时间：2022/7/15-2022/...
PhyloCSF.zip_coservation_lncRNA_phyloCSF分析_phylocsf值分析_rna保守性分析
09-23
总之，PhyloCSF是lncRNA保守性评估的重要工具，通过对多物种同源序列的分析，可以揭示lncRNA在进化上的保守性，为lncRNA的功能研究提供重要线索。"PhyloCSF.zip"中的"PhyloCSF.py"脚本则提供了实现这一分析的代码...
非线性非保守系统的稳定性判据
06-20
对于用Lienard方程描写的非线性自治电路,采用基波分析法,在适当端口施加正弦电压源uS,求得注入网络电流的基波分量IS1=Um(Gi+jBi).令基波输入导纳(Gi,Bi)=(0,0);如果求出有一组合理的实数解(ωS,Um)∈R2,说明网络存在...
基于海森矩阵的仿射谐波潮流保守性优化
01-15
针对传统仿射谐波潮流的计算结果存在较大保守性的问题，提出一种仿射谐波潮流的保守性优化方法。首先建立电网潮流方程，并采用牛顿-拉夫逊算法进行求解，在求解过程中采用纽曼级数将仿射矩阵的求逆运算转换为乘法...
【R语言】基于多模型的变量重要性图 (Variable Importance Plots)
amyniez的博客，欢迎交流讨论
08-22 1085
好久没有更新博客了，正好最近在帮老师做一个项目，里面涉及到了不同环境变量的重要性制图，所以在这里把我的理解进行分享，这应该是大家都可能遇到的问题。笔者水平有限，大家发现什么问题可以给我指出。变量重要度图(Variable importance plots)可以非常直观的展示各个变量在模型中的重要度，从而可以更好的理解和解释所建立的模型。
基于R语言的统计分析基础：数据结构
最新发布
m0_73500130的博客
08-25 592
R语言是一种用于统计分析和图形表示的编程语言和软件环境，它提供了多种数据结构以存储和操作数据。这些数据结构包括向量、矩阵、数组、数据框、列表、因子、Tibble、环境、公式、调用以及表达式。
基于R语言遥感随机森林建模与空间预测
weixin_48230888的博客
08-20 604
因此，遥感随机森林建模与空间预测的应用能够有效提升遥感数据分析的精度和可靠性，是许多研究者关注的热点。在R语言中，随机森林的实现与应用非常方便，R语言提供了多种包用于构建和优化随机森林模型。此外，R语言在数据可视化方面的优势使得用户能够直观地展示模型的结果和变量的重要性，进一步提高了分析的可解释性和应用价值。因此，R语言中的随机森林工具因其易用性、灵活性和强大的功能，成为遥感数据分析中不可或缺的工具。（2）R语言基础语法与数据结构，包括：程序包安装、加载、更新，数据读取与输出，ggplot2常规画图等。
R 语言学习教程，从入门到精通，R XML 文件使用（19）
qq_45746668的博客
08-20 316
XML 指的是可扩展标记语言（eXtensible Markup Language），XML 被设计用来传输和存储数据。 R 语言读写 XML 文件需要安装扩展包，我们可以在 R 到控制台输入以下命令来安装
广义时滞系统H∞控制：时滞相关与保守性分析
"该文章主要探讨了广义时滞系统的时滞相关型H∞控制问题，提出了新的Lyapunov-Krasovskii泛函和有界实引理，并设计了相应的控制器，减少了保守性。" 在控制系统理论中，时滞现象是实际系统中常见的一种非理想因素，...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值