view aurora_wgcna_trait.Rmd @ 7:ffbafe466107 draft

Uploaded

---
title: 'WGCNA: Module-Trait Association Report.'
output:
    html_document:
      number_sections: false
      toc: true
---

```{r setup, include=FALSE, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(error = TRUE, echo = FALSE)
```
```{r}
# Load the data from the previous step.
load(file=opt$r_data)
```
# Trait/Phenotype Data

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).

```{r}
# Load the trait data file.
trait_data = data.frame()
trait_data = read.csv(opt$trait_data, header = TRUE, row.names = opt$sname_col)
sample_names = rownames(gemt)
trait_rows = match(sample_names, rownames(trait_data))
trait_data = trait_data[trait_rows, ]

# Remove ignored columns.
ignore = strsplit(opt$ignore_cols, ',')
if (length(ignore[[1]]) > 0) {
  print(paste('You chose to ignore the following fields:', paste(ignore[[1]], collapse=", ")))
  trait_data = trait_data[, colnames(trait_data)[!(colnames(trait_data) %in% ignore[[1]])]]
}

# Change any categorical fields to 1 hot encoding as requested by the caller.
one_hot_cols = strsplit(opt$one_hot_cols, ',')
if (length(one_hot_cols[[1]]) > 0) {
  print(paste('You chose to 1-hot encode the following fields:', paste(one_hot_cols[[1]], collapse=", ")))

  # Perform the 1-hot encoding for specified fields.
  swap_cols = colnames(trait_data)[(colnames(trait_data) %in% one_hot_cols[[1]])]
  temp = as.data.frame(trait_data[, swap_cols])
  colnames(temp) = swap_cols
  temp = apply(temp, 2, make.names)
  dmy <- dummyVars(" ~ .", data = temp)
  encoded <- data.frame(predict(dmy, newdata = temp))
  encoded =  sapply(encoded, as.integer)

  # Make a new trait_data table with these new 1-hot fields.
  keep_cols = colnames(trait_data)[!(colnames(trait_data) %in% one_hot_cols[[1]])]
  keep = as.data.frame(trait_data[, keep_cols])
  colnames(keep) = keep_cols

  # Make a new trait_data object that has the columns to keep and the new 1-hot columns.
  trait_data = cbind(keep, encoded)
}

datatable(trait_data)
```
# Module-Condition Association.
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.

```{r fig.align='center', fig.width=8, fig.height=9}

# Determine the column types within the trait annotation data.
trait_types = sapply(trait_data, class)
trait_colors = data.frame(empty = rep(1:dim(trait_data)[1]))

# Set the colors for the quantitative data.
quantitative_fields = colnames(trait_data)[which(trait_types == "numeric")]
if (length(quantitative_fields) > 0) {
    qdata = as.data.frame(trait_data[,quantitative_fields])
    quantitative_colors = numbers2colors(qdata, signed = FALSE)
    colnames(quantitative_colors) = quantitative_fields
    trait_colors = cbind(trait_colors,quantitative_colors)
}

# Set the colors for the categorical data.
categorical_fields = colnames(trait_data)[which(trait_types == "factor")]
if (length(categorical_fields) > 0) {
    cdata = as.data.frame(trait_data[,categorical_fields])
    categorical_colors = labels2colors(cdata)
    colnames(categorical_colors) = categorical_fields
    trait_colors = cbind(trait_colors,categorical_colors)
}

# Set the colors for the ordinal data.
ordinal_fields = colnames(trait_data)[which(trait_types == "integer")]
if (length(ordinal_fields) > 0) {
    odata = as.data.frame(trait_data[,ordinal_fields])
    ordinal_colors = numbers2colors(odata, signed = FALSE)
    colnames(ordinal_colors) = ordinal_fields
    trait_colors = cbind(trait_colors, ordinal_colors)
}

trait_colors = subset(trait_colors, select=-c(empty))
trait_colors = trait_colors[,colnames(trait_data)]
options(repr.plot.width=15, repr.plot.height=10)
plotDendroAndColors(sampleTree, trait_colors,
                    groupLabels = names(trait_data),
                    main = "Sample Dendrogram and Annotation Heatmap",
                    cex.dendroLabels = 0.5)

```

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.

```{r fig.align='center', fig.width=15, fig.height=15}
MEs = orderMEs(MEs)
moduleTraitCor = cor(MEs, trait_data, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, n_samples);

# The WGCNA labeledHeatmap function is too overloaded with detail, we'll create a simpler plot.

plotData = melt(moduleTraitCor)
# We want to makes sure the order is the same as in the
# labeledHeatmap function (example above)
plotData$Var1 = factor(plotData$Var1, levels = rev(colnames(MEs)), ordered=TRUE)
# Now use ggplot2 to make a nicer image.
p <- ggplot(plotData, aes(Var2, Var1, fill=value)) +
  geom_tile() + xlab('Experimental Conditions') + ylab('WGCNA Modules') +
  scale_fill_gradient2(low = "#0072B2", high = "#D55E00",
                       mid = "white", midpoint = 0,
                       limit = c(-1,1), name="PCC") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1, size=15),
        axis.text.y = element_text(angle = 0, hjust=1, vjust=0.5, size=15),
        legend.text=element_text(size=15),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_blank())
print(p)
```

```{r}
output = cbind(moduleTraitCor, moduleTraitPvalue)
write.csv(output, file = opt$module_association_file, quote=FALSE, row.names=TRUE)
```
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.
```{r}
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(gemt, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), n_samples));
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

# Calculate the gene trait significance as a Pearson's R and p-value.
gts = as.data.frame(cor(gemt, trait_data, use = "p"));
gtsp = as.data.frame(corPvalueStudent(as.matrix(gts), n_samples));
colnames(gtsp) = c(paste("p", names(trait_data), sep="."))
colnames(gts) = c(paste("GS", names(trait_data), sep="."))

# Write out the gene information.
output = cbind(Module = module_labels, gts, gtsp)
write.csv(output, file = opt$gene_association_file, quote=FALSE, row.names=TRUE)

```
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.