comparison aurora_wgcna_trait.Rmd @ 7:ffbafe466107 draft

Uploaded

6:083886b3f3db

```
```{r}
# Load the data from the previous step.
load(file=opt$r_data)
```
# Sample Annotation Data

The table below shows the list of sample annotations you provided.  Sample annotations can include categorical data related to experimental conditions, or numerical trait or phenotype data for each sample.

```{r}
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, ]
datatable(trait_data)
```
# Module-Condition Association.
Now that we have sample annotations, 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 sample annotations 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)

# Set the colors for the quantitative data.
quantitative_fields = which(trait_types == "numeric")
quantitative_colors = numbers2colors(trait_data[,quantitative_fields], signed = FALSE)
colnames(quantitative_colors) = colnames(trait_data[,quantitative_fields])

# Set the colors for the categorical data.
categorical_fields = which(trait_types == "factor")
categorical_colors = labels2colors(trait_data[,categorical_fields])
colnames(categorical_colors) = colnames(trait_data[,categorical_fields])

# Set the colors for the ordinal data.
ordinal_fields = which(trait_types == "integer")
ordinal_colors = numbers2colors(trait_data[,ordinal_fields], signed = FALSE)
colnames(ordinal_colors) = colnames(trait_data[,ordinal_fields])

# Combine the colors in their original order as the trait_data data frame.
traitColors = cbind(quantitative_colors, categorical_colors, ordinal_colors)

plotDendroAndColors(sampleTree, traitColors,
                    groupLabels = names(trait_data),
                    main = "Sample Dendrogram and Annotation Heatmap",
                    cex.dendroLabels = 0.5)

```

```{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);

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(trait_data),
  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"))
```

```{r}
output = cbind(moduleTraitCor, moduleTraitPvalue)
write.csv(output, file = opt$module_association_file, quote=FALSE, row.names=TRUE)

# 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 sample annotations. 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 sample annotation features.

7:ffbafe466107

```
```{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)

```

```{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)

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