-
Notifications
You must be signed in to change notification settings - Fork 0
/
WGCNA3.R
141 lines (130 loc) · 5.63 KB
/
WGCNA3.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
getwd();
workingDir = "D:/thesis/proposal/codes/R_Code/wgcna/data";
setwd(workingDir);
# Load the WGCNA package
library(WGCNA)
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE);
# Load the expression and trait data saved in the first part
lnames = load(file = "..\\result\\cancer-healthy-01-dataInput.RData");
#The variable lnames contains the names of loaded variables.
lnames
# Load network data saved in the second part.
lnames = load(file = "..\\result\\cancer-healthy-02-dataInput.RData");
lnames
###############################################
####Quantifying module-trait associations######
###############################################
# Define numbers of genes and samples
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
sizeGrWindow(10,15)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
yLabels = names(MEs),
ySymbols = names(MEs),
colorLabels = FALSE,
colors = blueWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text = 0.5,
zlim = c(-1,1),
main = paste("Module-trait relationships"))
###############################################
####Gene Significance and Module Membership####
###############################################
# Define variable cancer containing the cancer column of datTrait
#cancer = as.data.frame(datTraits$cancer);
cancer = as.data.frame(datTraits$GastricCan);
names(cancer) = "cancer"
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datExpr, cancer, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(cancer), sep="");
names(GSPvalue) = paste("p.GS.", names(cancer), sep="");
###############################################
###########Intramodular analysis###############
###############################################
sizeGrWindow(15, 100);
par(mfrow = c(3,ceiling(length(modNames)/3)));
#### my code
impMir <- matrix(nrow = 0, ncol = 4)
colnames(impMir) <- c("module","miRNA","geneModuleMembership","geneTraitSignificance")
#### end of my code
for (module in modNames){
print(module)
column = match(module, modNames);
moduleGenes = moduleColors==module;
#### my code
miRNames = colnames(datExpr)[moduleGenes]
myGMM = abs(geneModuleMembership[moduleGenes, column])
myGTS = abs(geneTraitSignificance[moduleGenes, 1])
myIndex = myGMM > 0.7 & myGTS > 0.5
a <- matrix(nrow = 0, ncol = 4)
a <- cbind(module,miRNames[myIndex],myGMM[myIndex],myGTS[myIndex])
if(ncol(a)>1)
impMir <- rbind(impMir,a)
#### end of my code
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste(" Membership in", module),
ylab = "Gene sig for cancer",
main = paste(" \n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
}
write.table(impMir,"impMir.txt",sep = "\t", row.names = TRUE, col.names = TRUE)
###############################################
###########Summary output######################
###############################################
names(datExpr)
#names(datExpr)[moduleColors=="red"]
#annot = read.csv(file = "GeneAnnotation.csv");
# dim(annot)
# names(annot)
probes = names(datExpr)
# probes2annot = match(probes, annot$substanceBXH)
# The following is the number or probes without annotation:
# sum(is.na(probes2annot))
# Should return 0.
# Create the starting data frame
# geneInfo0 = data.frame(substanceBXH = probes,
# geneSymbol = annot$gene_symbol[probes2annot],
# LocusLinkID = annot$LocusLinkID[probes2annot],
# moduleColor = moduleColors,
# geneTraitSignificance,
# GSPvalue)
geneInfo0 = data.frame(substanceBXH = probes,
moduleColor = moduleColors,
geneTraitSignificance,
GSPvalue)
# Order modules by their significance for cancer
modOrder = order(-abs(cor(MEs, cancer, use = "p")));
# Add module membership information in the chosen order
for (mod in 1:ncol(geneModuleMembership))
{
oldNames = names(geneInfo0)
geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],
MMPvalue[, modOrder[mod]]);
names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
# Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance
geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.cancer));
geneInfo = geneInfo0[geneOrder, ]
write.csv(geneInfo, file = "..\\result\\geneInfo.csv")