|
3
|
1 ---
|
|
|
2 title: 'WGCNA: Module-Trait Association Report.'
|
|
|
3 output:
|
|
|
4 html_document:
|
|
|
5 number_sections: false
|
|
|
6 toc: true
|
|
|
7 ---
|
|
|
8
|
|
|
9 ```{r setup, include=FALSE, warning=FALSE, message=FALSE}
|
|
|
10 knitr::opts_chunk$set(error = TRUE, echo = FALSE)
|
|
|
11 ```
|
|
|
12 ```{r}
|
|
|
13 # Load the data from the previous step.
|
|
|
14 load(file=opt$r_data)
|
|
|
15 ```
|
|
7
|
16 # Trait/Phenotype Data
|
|
3
|
17
|
|
7
|
18 The table below shows the list of trait/phenotype data provided, but with any ignored columns removed and any categorical columns converted to a one-hot enconding (e.g. 0 when present 1 when not present).
|
|
3
|
19
|
|
|
20 ```{r}
|
|
7
|
21 # Load the trait data file.
|
|
3
|
22 trait_data = data.frame()
|
|
|
23 trait_data = read.csv(opt$trait_data, header = TRUE, row.names = opt$sname_col)
|
|
|
24 sample_names = rownames(gemt)
|
|
|
25 trait_rows = match(sample_names, rownames(trait_data))
|
|
|
26 trait_data = trait_data[trait_rows, ]
|
|
7
|
27
|
|
|
28 # Remove ignored columns.
|
|
|
29 ignore = strsplit(opt$ignore_cols, ',')
|
|
|
30 if (length(ignore[[1]]) > 0) {
|
|
|
31 print(paste('You chose to ignore the following fields:', paste(ignore[[1]], collapse=", ")))
|
|
|
32 trait_data = trait_data[, colnames(trait_data)[!(colnames(trait_data) %in% ignore[[1]])]]
|
|
|
33 }
|
|
|
34
|
|
|
35 # Change any categorical fields to 1 hot encoding as requested by the caller.
|
|
|
36 one_hot_cols = strsplit(opt$one_hot_cols, ',')
|
|
|
37 if (length(one_hot_cols[[1]]) > 0) {
|
|
|
38 print(paste('You chose to 1-hot encode the following fields:', paste(one_hot_cols[[1]], collapse=", ")))
|
|
|
39
|
|
|
40 # Perform the 1-hot encoding for specified fields.
|
|
|
41 swap_cols = colnames(trait_data)[(colnames(trait_data) %in% one_hot_cols[[1]])]
|
|
|
42 temp = as.data.frame(trait_data[, swap_cols])
|
|
|
43 colnames(temp) = swap_cols
|
|
|
44 temp = apply(temp, 2, make.names)
|
|
|
45 dmy <- dummyVars(" ~ .", data = temp)
|
|
|
46 encoded <- data.frame(predict(dmy, newdata = temp))
|
|
|
47 encoded = sapply(encoded, as.integer)
|
|
|
48
|
|
|
49 # Make a new trait_data table with these new 1-hot fields.
|
|
|
50 keep_cols = colnames(trait_data)[!(colnames(trait_data) %in% one_hot_cols[[1]])]
|
|
|
51 keep = as.data.frame(trait_data[, keep_cols])
|
|
|
52 colnames(keep) = keep_cols
|
|
|
53
|
|
|
54 # Make a new trait_data object that has the columns to keep and the new 1-hot columns.
|
|
|
55 trait_data = cbind(keep, encoded)
|
|
|
56 }
|
|
|
57
|
|
3
|
58 datatable(trait_data)
|
|
|
59 ```
|
|
|
60 # Module-Condition Association.
|
|
7
|
61 Now that we have trait/phenotype data, we can explore if any of the network modules are asociated with these features. First, is an empirical exploration by viewing again the sample dendrogram but with traits added and colored by category or numerical intensity, as appropriate. If groups of samples with similar expression also share similar annotations then the same colors will appear "in blocks" under the clustered samples. This view does not indicate associations but can help visualize when some modules might be associated.
|
|
3
|
62
|
|
|
63 ```{r fig.align='center', fig.width=8, fig.height=9}
|
|
|
64
|
|
|
65 # Determine the column types within the trait annotation data.
|
|
|
66 trait_types = sapply(trait_data, class)
|
|
7
|
67 trait_colors = data.frame(empty = rep(1:dim(trait_data)[1]))
|
|
3
|
68
|
|
|
69 # Set the colors for the quantitative data.
|
|
7
|
70 quantitative_fields = colnames(trait_data)[which(trait_types == "numeric")]
|
|
|
71 if (length(quantitative_fields) > 0) {
|
|
|
72 qdata = as.data.frame(trait_data[,quantitative_fields])
|
|
|
73 quantitative_colors = numbers2colors(qdata, signed = FALSE)
|
|
|
74 colnames(quantitative_colors) = quantitative_fields
|
|
|
75 trait_colors = cbind(trait_colors,quantitative_colors)
|
|
|
76 }
|
|
3
|
77
|
|
|
78 # Set the colors for the categorical data.
|
|
7
|
79 categorical_fields = colnames(trait_data)[which(trait_types == "factor")]
|
|
|
80 if (length(categorical_fields) > 0) {
|
|
|
81 cdata = as.data.frame(trait_data[,categorical_fields])
|
|
|
82 categorical_colors = labels2colors(cdata)
|
|
|
83 colnames(categorical_colors) = categorical_fields
|
|
|
84 trait_colors = cbind(trait_colors,categorical_colors)
|
|
|
85 }
|
|
3
|
86
|
|
|
87 # Set the colors for the ordinal data.
|
|
7
|
88 ordinal_fields = colnames(trait_data)[which(trait_types == "integer")]
|
|
|
89 if (length(ordinal_fields) > 0) {
|
|
|
90 odata = as.data.frame(trait_data[,ordinal_fields])
|
|
|
91 ordinal_colors = numbers2colors(odata, signed = FALSE)
|
|
|
92 colnames(ordinal_colors) = ordinal_fields
|
|
|
93 trait_colors = cbind(trait_colors, ordinal_colors)
|
|
|
94 }
|
|
3
|
95
|
|
7
|
96 trait_colors = subset(trait_colors, select=-c(empty))
|
|
|
97 trait_colors = trait_colors[,colnames(trait_data)]
|
|
|
98 options(repr.plot.width=15, repr.plot.height=10)
|
|
|
99 plotDendroAndColors(sampleTree, trait_colors,
|
|
3
|
100 groupLabels = names(trait_data),
|
|
|
101 main = "Sample Dendrogram and Annotation Heatmap",
|
|
|
102 cex.dendroLabels = 0.5)
|
|
|
103
|
|
|
104 ```
|
|
|
105
|
|
|
106 To statistically identify the associations, correlation tests are performed of the eigengenes of each module with the annotation data. The following heatmap shows the results between each annotation feature and each module. Modules with a signficant positive assocation have a correlation value near 1. Modules with a significant negative association have a correlation value near -1. Modules with no correlation have a value near 0.
|
|
|
107
|
|
|
108 ```{r fig.align='center', fig.width=15, fig.height=15}
|
|
|
109 MEs = orderMEs(MEs)
|
|
|
110 moduleTraitCor = cor(MEs, trait_data, use = "p");
|
|
|
111 moduleTraitPvalue = corPvalueStudent(moduleTraitCor, n_samples);
|
|
|
112
|
|
7
|
113 # The WGCNA labeledHeatmap function is too overloaded with detail, we'll create a simpler plot.
|
|
3
|
114
|
|
7
|
115 plotData = melt(moduleTraitCor)
|
|
|
116 # We want to makes sure the order is the same as in the
|
|
|
117 # labeledHeatmap function (example above)
|
|
|
118 plotData$Var1 = factor(plotData$Var1, levels = rev(colnames(MEs)), ordered=TRUE)
|
|
|
119 # Now use ggplot2 to make a nicer image.
|
|
|
120 p <- ggplot(plotData, aes(Var2, Var1, fill=value)) +
|
|
|
121 geom_tile() + xlab('Experimental Conditions') + ylab('WGCNA Modules') +
|
|
|
122 scale_fill_gradient2(low = "#0072B2", high = "#D55E00",
|
|
|
123 mid = "white", midpoint = 0,
|
|
|
124 limit = c(-1,1), name="PCC") +
|
|
|
125 theme_bw() +
|
|
|
126 theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1, size=15),
|
|
|
127 axis.text.y = element_text(angle = 0, hjust=1, vjust=0.5, size=15),
|
|
|
128 legend.text=element_text(size=15),
|
|
|
129 panel.border = element_blank(),
|
|
|
130 panel.grid.major = element_blank(),
|
|
|
131 panel.grid.minor = element_blank(),
|
|
|
132 axis.line = element_blank())
|
|
|
133 print(p)
|
|
3
|
134 ```
|
|
|
135
|
|
|
136 ```{r}
|
|
|
137 output = cbind(moduleTraitCor, moduleTraitPvalue)
|
|
|
138 write.csv(output, file = opt$module_association_file, quote=FALSE, row.names=TRUE)
|
|
|
139 ```
|
|
|
140 A file has been generated named `module_association.csv` which conatins the list of modules, and their correlation values as well as p-values indicating the strength of the associations.
|
|
|
141 ```{r}
|
|
|
142 # names (colors) of the modules
|
|
|
143 modNames = substring(names(MEs), 3)
|
|
|
144 geneModuleMembership = as.data.frame(cor(gemt, MEs, use = "p"));
|
|
|
145 MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), n_samples));
|
|
|
146 names(geneModuleMembership) = paste("MM", modNames, sep="");
|
|
|
147 names(MMPvalue) = paste("p.MM", modNames, sep="");
|
|
|
148
|
|
|
149 # Calculate the gene trait significance as a Pearson's R and p-value.
|
|
|
150 gts = as.data.frame(cor(gemt, trait_data, use = "p"));
|
|
|
151 gtsp = as.data.frame(corPvalueStudent(as.matrix(gts), n_samples));
|
|
|
152 colnames(gtsp) = c(paste("p", names(trait_data), sep="."))
|
|
|
153 colnames(gts) = c(paste("GS", names(trait_data), sep="."))
|
|
|
154
|
|
|
155 # Write out the gene information.
|
|
|
156 output = cbind(Module = module_labels, gts, gtsp)
|
|
|
157 write.csv(output, file = opt$gene_association_file, quote=FALSE, row.names=TRUE)
|
|
|
158
|
|
|
159 ```
|
|
7
|
160 Genes themselves can also have assocation with traits. This is calculated via a traditional correlation test as well. Another file has been generated named `gene_association.csv` which provides the list of genes, the modules they belong to and the assocaition of each gene to the trait features.
|
