# HG changeset patch
# User mingchen0919
# Date 1502210150 14400
# Node ID 4275479ada3a44d1f94cc23fe10c0b29c909219e
planemo upload for repository https://github.com/statonlab/docker-GRReport/tree/master/my_tools/rmarkdown_wgcna commit d91f269e8bc09a488ed2e005122bbb4a521f44a0-dirty
diff -r 000000000000 -r 4275479ada3a wgcna_construct_network.Rmd
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_construct_network.Rmd Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,178 @@
+---
+title: 'WGCNA: construct network'
+output:
+ html_document:
+ number_sections: true
+ toc: true
+ theme: cosmo
+ highlight: tango
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+ echo = ECHO
+)
+```
+
+# Import workspace
+
+This step imports workspace from the **WGCNA: preprocessing** step.
+
+```{r}
+fcp = file.copy("PREPROCESSING_WORKSPACE", "deseq.RData")
+load("deseq.RData")
+```
+
+
+# Processing outliers {.tabset}
+
+## Before removing outliers
+
+```{r}
+plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
+ cex.axis = 1, cex.main = 1, cex = 0.5)
+if(!is.na(HEIGHT_CUT)) {
+ # plot a line to show the cut
+ abline(h = HEIGHT_CUT, col = "red")
+ # determine cluster under the line
+ clust = cutreeStatic(sampleTree, cutHeight = HEIGHT_CUT, minSize = 10)
+ keepSamples = (clust==1)
+ expression_data = expression_data[keepSamples, ]
+}
+```
+
+## After removing outliers
+
+```{r}
+sampleTree = hclust(dist(expression_data), method = "average");
+plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="",
+ cex.axis = 1, cex.main = 1, cex = 0.5)
+```
+
+
+# Trait data {.tabeset}
+
+If trait data is provided, the first 100 rows from the data will be displayed here. A plot consisting of sample cluster dendrogram and trait heatmap will also be gerenated.
+
+## Trait data table
+
+```{r}
+trait_data = data.frame()
+if ("TRAIT_DATA" != 'None') {
+ trait_data = read.csv("TRAIT_DATA", header = TRUE, row.names = 1)
+ # form a data frame analogous to expression data that will hold the traits.
+ sample_names = rownames(expression_data)
+ trait_rows = match(sample_names, rownames(trait_data))
+ trait_data = trait_data[trait_rows, ]
+ datatable(head(trait_data, 100), style="bootstrap", filter = 'top',
+ class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+}
+```
+
+## Dendrogram and heatmap
+
+```{r fig.align='center', fig.width=8, fig.height=9}
+if (nrow(trait_data) != 0) {
+ traitColors = numbers2colors(trait_data, signed = FALSE)
+ plotDendroAndColors(sampleTree, traitColors,
+ groupLabels = names(trait_data),
+ main = "Sample dendrogram and trait heatmap",
+ cex.dendroLabels = 0.5)
+}
+```
+
+
+# The thresholding power
+
+```{r}
+powers = c(1:10, seq(12, 20, 2))
+soft_threshold = pickSoftThreshold(expression_data, powerVector = powers, verbose = 5)
+```
+
+```{r fig.align='center'}
+par(mfrow=c(1,2))
+plot(soft_threshold$fitIndices[,1], -sign(soft_threshold$fitIndices[,3])*soft_threshold$fitIndices[,2],
+ xlab="Soft Threshold (power)",
+ ylab="Scale Free Topology Model Fit,signed R^2",type="n",
+ main = paste("Scale independence"),
+ cex.lab = 0.5);
+text(soft_threshold$fitIndices[,1], -sign(soft_threshold$fitIndices[,3])*soft_threshold$fitIndices[,2],
+ labels=powers,cex=0.5,col="red");
+
+# calculate soft threshold power
+y = -sign(soft_threshold$fitIndices[,3])*soft_threshold$fitIndices[,2]
+r2_cutoff = 0.9
+for(i in 1:length(powers)) {
+ if(y[i] > r2_cutoff) {
+ soft_threshold_power = soft_threshold$fitIndices[,1][i]
+ r2_cutoff_new = y[i]
+ break
+ }
+ soft_threshold_power = soft_threshold$fitIndices[,1][length(powers)]
+}
+abline(h=r2_cutoff, col="red")
+abline(v=soft_threshold_power, col="blue")
+text(soft_threshold_power+1, r2_cutoff-0.1,
+ paste0('R^2 cutoff = ', round(r2_cutoff_new,2)),
+ cex = 0.5, col = "red")
+
+plot(soft_threshold$fitIndices[,1], soft_threshold$fitIndices[,5],
+ xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
+ main = paste("Mean connectivity"),
+ cex.lab = 0.5)
+text(soft_threshold$fitIndices[,1], soft_threshold$fitIndices[,5], labels=powers, cex=0.5,col="red")
+par(mfrow=c(1,1))
+```
+
+
+# Construct network
+
+The gene network is constructed based on **soft threshold power = `r soft_threshold_power`**
+
+```{r}
+gene_network = blockwiseModules(expression_data, power = soft_threshold_power,
+ TOMType = "unsigned", minModuleSize = 30,
+ reassignThreshold = 0, mergeCutHeight = 0.25,
+ numericLabels = TRUE, pamRespectsDendro = FALSE,
+ verbose = 3)
+```
+
+
+# Gene modules {.tabset}
+
+## Idenfity gene modules
+
+```{r}
+modules = table(gene_network$colors)
+n_modules = length(modules) - 1
+module_size_upper = modules[2]
+module_size_lower = modules[length(modules)]
+
+module_table = data.frame(model_label = c(0, 1:n_modules),
+ gene_size = as.vector(modules))
+datatable(t(module_table))
+```
+
+The results above indicates that there are **`r n_modules` gene modules**, labeled 1 through `r length(n_modules)` in order of descending size. The largest module has **`r module_size_upper` genes**, and the smallest module has **`r module_size_lower` genes**. The label 0 is reserved for genes outside of all modules.
+
+
+## Dendrogram and module plot
+
+```{r}
+# Convert labels to colors for plotting
+module_colors = labels2colors(gene_network$colors)
+# Plot the dendrogram and the module colors underneath
+plotDendroAndColors(gene_network$dendrograms[[1]], module_colors[gene_network$blockGenes[[1]]],
+ "Module colors",
+ dendroLabels = FALSE, hang = 0.03,
+ addGuide = TRUE, guideHang = 0.05)
+```
+
+
+```{r echo=FALSE}
+# save workspace
+rm("opt")
+save(list=ls(all.names = TRUE), file='CONSTRUCT_NETWORK_WORKSPACE')
+```
+
+
diff -r 000000000000 -r 4275479ada3a wgcna_construct_network.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_construct_network.xml Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,105 @@
+
+
+ r-getopt
+ r-rmarkdown
+ r-plyr
+ r-highcharter
+ r-dt
+ r-htmltools
+ r-wgcna
+
+
+ Construct gene network.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ @article{langfelder2008wgcna,
+ title={WGCNA: an R package for weighted correlation network analysis},
+ author={Langfelder, Peter and Horvath, Steve},
+ journal={BMC bioinformatics},
+ volume={9},
+ number={1},
+ pages={559},
+ year={2008},
+ publisher={BioMed Central}
+ }
+
+
+ @article{allaire2016rmarkdown,
+ title={rmarkdown: Dynamic Documents for R, 2016},
+ author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff and Wickham, Hadley and Atkins, Aron and Hyndman, Rob},
+ journal={R package version 0.9},
+ volume={6},
+ year={2016}
+ }
+
+
+ @book{xie2015dynamic,
+ title={Dynamic Documents with R and knitr},
+ author={Xie, Yihui},
+ volume={29},
+ year={2015},
+ publisher={CRC Press}
+ }
+
+
+
\ No newline at end of file
diff -r 000000000000 -r 4275479ada3a wgcna_construct_network_render.R
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_construct_network_render.R Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,112 @@
+##======= Handle arguments from command line ========
+# setup R error handline to go to stderr
+options(show.error.messages=FALSE,
+ error=function(){
+ cat(geterrmessage(), file=stderr())
+ quit("no", 1, F)
+ })
+
+# we need that to not crash galaxy with an UTF8 error on German LC settings.
+loc = Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+
+# suppress warning
+options(warn = -1)
+
+options(stringsAsFactors=FALSE, useFancyQuotes=FALSE)
+args = commandArgs(trailingOnly=TRUE)
+
+suppressPackageStartupMessages({
+ library(getopt)
+ library(tools)
+})
+
+# column 1: the long flag name
+# column 2: the short flag alias. A SINGLE character string
+# column 3: argument mask
+# 0: no argument
+# 1: argument required
+# 2: argument is optional
+# column 4: date type to which the flag's argument shall be cast.
+# possible values: logical, integer, double, complex, character.
+spec_list=list()
+
+##------- 1. input data ---------------------
+spec_list$ECHO = c('echo', 'e', '1', 'character')
+spec_list$PREPROCESSING_WORKSPACE = c('preprocessing_workspace', 'w', '1', 'character')
+spec_list$HEIGHT_CUT = c('height_cut', 'h', '2', 'double')
+spec_list$TRAIT_DATA = c('trait_data', 't', '2', 'character')
+
+
+##--------2. output report and report site directory --------------
+spec_list$OUTPUT_HTML = c('wgcna_construct_network_html', 'o', '1', 'character')
+spec_list$OUTPUT_DIR = c('wgcna_construct_network_dir', 'd', '1', 'character')
+spec_list$CONSTRUCT_NETWORK_WORKSPACE = c('construct_network_workspace', 'W', '1', 'character')
+
+
+##--------3. Rmd templates in the tool directory ----------
+
+spec_list$WGCNA_PREPROCESSING_RMD = c('wgcna_construct_network_rmd', 'M', '1', 'character')
+
+
+
+##------------------------------------------------------------------
+
+spec = t(as.data.frame(spec_list))
+opt = getopt(spec)
+# arguments are accessed by long flag name (the first column in the spec matrix)
+# NOT by element name in the spec_list
+# example: opt$help, opt$expression_file
+##====== End of arguments handling ==========
+
+#------ Load libraries ---------
+library(rmarkdown)
+library(WGCNA)
+library(DT)
+library(htmltools)
+library(ggplot2)
+
+
+#----- 1. create the report directory ------------------------
+system(paste0('mkdir -p ', opt$wgcna_construct_network_dir))
+
+
+#----- 2. generate Rmd files with Rmd templates --------------
+# a. templates without placeholder variables:
+# copy templates from tool directory to the working directory.
+# b. templates with placeholder variables:
+# substitute variables with user input values and place them in the working directory.
+
+
+#----- 01 wgcna_construct_network.Rmd -----------------------
+readLines(opt$wgcna_construct_network_rmd) %>%
+ (function(x) {
+ gsub('ECHO', opt$echo, x)
+ }) %>%
+ (function(x) {
+ gsub('PREPROCESSING_WORKSPACE', opt$preprocessing_workspace, x)
+ }) %>%
+ (function(x) {
+ gsub('HEIGHT_CUT', opt$height_cut, x)
+ }) %>%
+ (function(x) {
+ gsub('TRAIT_DATA', opt$trait_data, x)
+ }) %>%
+ (function(x) {
+ gsub('OUTPUT_DIR', opt$wgcna_construct_network_dir, x)
+ }) %>%
+ (function(x) {
+ gsub('CONSTRUCT_NETWORK_WORKSPACE', opt$construct_network_workspace, x)
+ }) %>%
+ (function(x) {
+ fileConn = file('wgcna_construct_network.Rmd')
+ writeLines(x, con=fileConn)
+ close(fileConn)
+ })
+
+
+#------ 3. render all Rmd files --------
+render('wgcna_construct_network.Rmd', output_file = opt$wgcna_construct_network_html)
+
+#-------4. manipulate outputs -----------------------------
+
+
diff -r 000000000000 -r 4275479ada3a wgcna_eigengene_visualization.Rmd
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_eigengene_visualization.Rmd Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,121 @@
+---
+title: 'WGCNA: eigengene visualization'
+output:
+ html_document:
+ number_sections: true
+ toc: true
+ theme: cosmo
+ highlight: tango
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+ echo = ECHO
+)
+```
+
+# Import workspace
+
+This step imports workspace from the **WGCNA: construct network** step.
+
+```{r}
+fcp = file.copy("CONSTRUCT_NETWORK_WORKSPACE", "deseq.RData")
+load("deseq.RData")
+```
+
+
+# Gene modules {.tabset}
+
+```{r}
+if(!is.na(SOFT_THRESHOLD_POWER)) soft_threshold_power = SOFT_THRESHOLD_POWER
+```
+
+## Identify gene modules
+
+The gene network is constructed based on **soft threshold power = `r soft_threshold_power`**
+
+```{r}
+gene_network = blockwiseModules(expression_data, power = soft_threshold_power,
+ TOMType = "unsigned", minModuleSize = 30,
+ reassignThreshold = 0, mergeCutHeight = 0.25,
+ numericLabels = TRUE, pamRespectsDendro = FALSE,
+ verbose = 3)
+```
+
+
+```{r}
+modules = table(gene_network$colors)
+n_modules = length(modules) - 1
+module_size_upper = modules[2]
+module_size_lower = modules[length(modules)]
+
+module_table = data.frame(model_label = c(0, 1:n_modules),
+ gene_size = as.vector(modules))
+datatable(t(module_table))
+```
+
+The results above indicates that there are **`r n_modules` gene modules**, labeled 1 through `r length(n_modules)` in order of descending size. The largest module has **`r module_size_upper` genes**, and the smallest module has **`r module_size_lower` genes**. The label 0 is reserved for genes outside of all modules.
+
+
+## Dendrogram and module plot
+
+```{r}
+# Convert labels to colors for plotting
+module_colors = labels2colors(gene_network$colors)
+# Plot the dendrogram and the module colors underneath
+plotDendroAndColors(gene_network$dendrograms[[1]], module_colors[gene_network$blockGenes[[1]]],
+ "Module colors",
+ dendroLabels = FALSE, hang = 0.03,
+ addGuide = TRUE, guideHang = 0.05)
+```
+
+
+# Gene module correlation
+
+We can calculate eigengenes and use them as representative profiles to quantify similarity of found gene modules.
+
+```{r}
+n_genes = ncol(expression_data)
+n_samples = nrow(expression_data)
+```
+
+```{r}
+diss_tom = 1-TOMsimilarityFromExpr(expression_data, power = soft_threshold_power)
+set.seed(123)
+select_genes = sample(n_genes, size = PLOT_GENES)
+select_diss_tom = diss_tom[select_genes, select_genes]
+
+# calculate gene tree on selected genes
+select_gene_tree = hclust(as.dist(select_diss_tom), method = 'average')
+select_module_colors = module_colors[select_genes]
+
+# transform diss_tom with a power to make moderately strong connections more visiable in the heatmap.
+plot_diss_tom = select_diss_tom^7
+# set diagonal to NA for a nicer plot
+diag(plot_diss_tom) = NA
+```
+
+
+```{r fig.align='center'}
+TOMplot(plot_diss_tom, select_gene_tree, select_module_colors, main = "Network heatmap")
+```
+
+
+# Eigengene visualization {.tabset}
+
+## Eigengene dendrogram
+
+```{r fig.align='center'}
+module_eigengenes = moduleEigengenes(expression_data, module_colors)$eigengenes
+plotEigengeneNetworks(module_eigengenes, "Eigengene dendrogram",
+ plotHeatmaps = FALSE)
+```
+
+## Eigengene adjacency heatmap
+
+```{r fig.align='center'}
+plotEigengeneNetworks(module_eigengenes, "Eigengene adjacency heatmap",
+ marHeatmap = c(2, 3, 2, 2),
+ plotDendrograms = FALSE, xLabelsAngle = 90)
+```
+
diff -r 000000000000 -r 4275479ada3a wgcna_eigengene_visualization.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_eigengene_visualization.xml Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,100 @@
+
+
+ r-getopt
+ r-rmarkdown
+ r-plyr
+ r-highcharter
+ r-dt
+ r-htmltools
+ r-wgcna
+
+
+ Eigengene visualization.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ @article{langfelder2008wgcna,
+ title={WGCNA: an R package for weighted correlation network analysis},
+ author={Langfelder, Peter and Horvath, Steve},
+ journal={BMC bioinformatics},
+ volume={9},
+ number={1},
+ pages={559},
+ year={2008},
+ publisher={BioMed Central}
+ }
+
+
+ @article{allaire2016rmarkdown,
+ title={rmarkdown: Dynamic Documents for R, 2016},
+ author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff and Wickham, Hadley and Atkins, Aron and Hyndman, Rob},
+ journal={R package version 0.9},
+ volume={6},
+ year={2016}
+ }
+
+
+ @book{xie2015dynamic,
+ title={Dynamic Documents with R and knitr},
+ author={Xie, Yihui},
+ volume={29},
+ year={2015},
+ publisher={CRC Press}
+ }
+
+
+
\ No newline at end of file
diff -r 000000000000 -r 4275479ada3a wgcna_eigengene_visualization_render.R
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_eigengene_visualization_render.R Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,109 @@
+##======= Handle arguments from command line ========
+# setup R error handline to go to stderr
+options(show.error.messages=FALSE,
+ error=function(){
+ cat(geterrmessage(), file=stderr())
+ quit("no", 1, F)
+ })
+
+# we need that to not crash galaxy with an UTF8 error on German LC settings.
+loc = Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+
+# suppress warning
+options(warn = -1)
+
+options(stringsAsFactors=FALSE, useFancyQuotes=FALSE)
+args = commandArgs(trailingOnly=TRUE)
+
+suppressPackageStartupMessages({
+ library(getopt)
+ library(tools)
+})
+
+# column 1: the long flag name
+# column 2: the short flag alias. A SINGLE character string
+# column 3: argument mask
+# 0: no argument
+# 1: argument required
+# 2: argument is optional
+# column 4: date type to which the flag's argument shall be cast.
+# possible values: logical, integer, double, complex, character.
+spec_list=list()
+
+##------- 1. input data ---------------------
+spec_list$ECHO = c('echo', 'e', '1', 'character')
+spec_list$CONSTRUCT_NETWORK_WORKSPACE = c('construct_network_workspace', 'w', '1', 'character')
+spec_list$SOFT_THRESHOLD_POWER = c('soft_threshold_power', 'p', '2', 'double')
+spec_list$PLOT_GENES = c('plot_genes', 'n', '1', 'integer')
+
+
+##--------2. output report and report site directory --------------
+spec_list$OUTPUT_HTML = c('wgcna_eigengene_visualization_html', 'o', '1', 'character')
+spec_list$OUTPUT_DIR = c('wgcna_eigengene_visualization_dir', 'd', '1', 'character')
+
+
+
+##--------3. Rmd templates in the tool directory ----------
+
+spec_list$WGCNA_EIGENGENE_VISUALIZATION_RMD = c('wgcna_eigengene_visualization_rmd', 'M', '1', 'character')
+
+
+
+##------------------------------------------------------------------
+
+spec = t(as.data.frame(spec_list))
+opt = getopt(spec)
+# arguments are accessed by long flag name (the first column in the spec matrix)
+# NOT by element name in the spec_list
+# example: opt$help, opt$expression_file
+##====== End of arguments handling ==========
+
+#------ Load libraries ---------
+library(rmarkdown)
+library(WGCNA)
+library(DT)
+library(htmltools)
+library(ggplot2)
+
+
+#----- 1. create the report directory ------------------------
+system(paste0('mkdir -p ', opt$wgcna_eigengene_visualization_dir))
+
+
+#----- 2. generate Rmd files with Rmd templates --------------
+# a. templates without placeholder variables:
+# copy templates from tool directory to the working directory.
+# b. templates with placeholder variables:
+# substitute variables with user input values and place them in the working directory.
+
+
+#----- 01 wgcna_eigengene_visualization.Rmd -----------------------
+readLines(opt$wgcna_eigengene_visualization_rmd) %>%
+ (function(x) {
+ gsub('ECHO', opt$echo, x)
+ }) %>%
+ (function(x) {
+ gsub('CONSTRUCT_NETWORK_WORKSPACE', opt$construct_network_workspace, x)
+ }) %>%
+ (function(x) {
+ gsub('SOFT_THRESHOLD_POWER', opt$soft_threshold_power, x)
+ }) %>%
+ (function(x) {
+ gsub('PLOT_GENES', opt$plot_genes, x)
+ }) %>%
+ (function(x) {
+ gsub('OUTPUT_DIR', opt$wgcna_eigengene_visualization_dir, x)
+ }) %>%
+ (function(x) {
+ fileConn = file('wgcna_eigengene_visualization.Rmd')
+ writeLines(x, con=fileConn)
+ close(fileConn)
+ })
+
+
+#------ 3. render all Rmd files --------
+render('wgcna_eigengene_visualization.Rmd', output_file = opt$wgcna_eigengene_visualization_html)
+
+#-------4. manipulate outputs -----------------------------
+
+
diff -r 000000000000 -r 4275479ada3a wgcna_preprocessing.Rmd
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_preprocessing.Rmd Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,76 @@
+---
+title: 'WGCNA: data preprocessing'
+output:
+ html_document:
+ number_sections: true
+ toc: true
+ theme: cosmo
+ highlight: tango
+---
+
+```{r setup, include=FALSE, warning=FALSE, message=FALSE}
+knitr::opts_chunk$set(
+ echo = ECHO
+)
+```
+
+```{r}
+str(opt)
+```
+
+# Import data
+
+Each row represents a gene and each column represents a sample.
+
+```{r}
+expression_data = read.csv('EXPRESSION_DATA', header = TRUE, row.names = 1)
+```
+
+Display the first 100 genes.
+
+```{r}
+datatable(head(expression_data, 100), style="bootstrap", filter = 'top',
+ class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
+```
+
+Transpose expression data matrix so that each row represents a sample and each column represents a gene.
+
+```{r}
+expression_data = as.data.frame(t(expression_data))
+```
+
+# Checking data
+
+Checking data for excessive missing values and identification of outlier microarray samples.
+
+```{r}
+gsg = goodSamplesGenes(expression_data, verbose = 3)
+if (!gsg$allOK) {
+ # Optionally, print the gene and sample names that were removed:
+ if (sum(!gsg$goodGenes)>0)
+ printFlush(paste("Removing genes:", paste(names(expression_data)[!gsg$goodGenes], collapse = ", ")));
+ if (sum(!gsg$goodSamples)>0)
+ printFlush(paste("Removing samples:", paste(rownames(expression_data)[!gsg$goodSamples], collapse = ", ")));
+ # Remove the offending genes and samples from the data:
+ expression_data = expression_data[gsg$goodSamples, gsg$goodGenes]
+} else {
+ print('all genes are OK!')
+}
+```
+
+# Clustering samples
+
+If there are any outliers, choose a height cut that will remove the offending sample. Remember this number since you will need this number in further analysis.
+
+```{r fig.align='center'}
+sampleTree = hclust(dist(expression_data), method = "average");
+plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="",
+ cex.axis = 1, cex.main = 1, cex = 0.5)
+```
+
+
+```{r echo=FALSE}
+rm("opt")
+save(list=ls(all.names = TRUE), file='PREPROCESSING_WORKSPACE')
+```
+
diff -r 000000000000 -r 4275479ada3a wgcna_preprocessing.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_preprocessing.xml Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,96 @@
+
+
+ r-getopt
+ r-rmarkdown
+ r-plyr
+ r-highcharter
+ r-dt
+ r-htmltools
+ r-wgcna
+
+
+ Data clearning and preprocessing.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ @article{langfelder2008wgcna,
+ title={WGCNA: an R package for weighted correlation network analysis},
+ author={Langfelder, Peter and Horvath, Steve},
+ journal={BMC bioinformatics},
+ volume={9},
+ number={1},
+ pages={559},
+ year={2008},
+ publisher={BioMed Central}
+ }
+
+
+ @article{allaire2016rmarkdown,
+ title={rmarkdown: Dynamic Documents for R, 2016},
+ author={Allaire, J and Cheng, Joe and Xie, Yihui and McPherson, Jonathan and Chang, Winston and Allen, Jeff and Wickham, Hadley and Atkins, Aron and Hyndman, Rob},
+ journal={R package version 0.9},
+ volume={6},
+ year={2016}
+ }
+
+
+ @book{xie2015dynamic,
+ title={Dynamic Documents with R and knitr},
+ author={Xie, Yihui},
+ volume={29},
+ year={2015},
+ publisher={CRC Press}
+ }
+
+
+
\ No newline at end of file
diff -r 000000000000 -r 4275479ada3a wgcna_preprocessing_render.R
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/wgcna_preprocessing_render.R Tue Aug 08 12:35:50 2017 -0400
@@ -0,0 +1,102 @@
+##======= Handle arguments from command line ========
+# setup R error handline to go to stderr
+options(show.error.messages=FALSE,
+ error=function(){
+ cat(geterrmessage(), file=stderr())
+ quit("no", 1, F)
+ })
+
+# we need that to not crash galaxy with an UTF8 error on German LC settings.
+loc = Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
+
+# suppress warning
+options(warn = -1)
+
+options(stringsAsFactors=FALSE, useFancyQuotes=FALSE)
+args = commandArgs(trailingOnly=TRUE)
+
+suppressPackageStartupMessages({
+ library(getopt)
+ library(tools)
+})
+
+# column 1: the long flag name
+# column 2: the short flag alias. A SINGLE character string
+# column 3: argument mask
+# 0: no argument
+# 1: argument required
+# 2: argument is optional
+# column 4: date type to which the flag's argument shall be cast.
+# possible values: logical, integer, double, complex, character.
+spec_list=list()
+
+##------- 1. input data ---------------------
+spec_list$ECHO = c('echo', 'e', '1', 'character')
+spec_list$EXPRESSION_DATA = c('expression_data', 'E', '1', 'character')
+
+
+##--------2. output report and report site directory --------------
+spec_list$OUTPUT_HTML = c('wgcna_preprocessing_html', 'o', '1', 'character')
+spec_list$OUTPUT_DIR = c('wgcna_preprocessing_dir', 'd', '1', 'character')
+spec_list$PREPROCESSING_WORKSPACE = c('preprocessing_workspace', 'w', '1', 'character')
+
+##--------3. Rmd templates sitting in the tool directory ----------
+
+spec_list$WGCNA_PREPROCESSING_RMD = c('wgcna_preprocessing_rmd', 'D', '1', 'character')
+
+
+
+##------------------------------------------------------------------
+
+spec = t(as.data.frame(spec_list))
+opt = getopt(spec)
+# arguments are accessed by long flag name (the first column in the spec matrix)
+# NOT by element name in the spec_list
+# example: opt$help, opt$expression_file
+##====== End of arguments handling ==========
+
+#------ Load libraries ---------
+library(rmarkdown)
+library(WGCNA)
+library(DT)
+library(htmltools)
+
+
+#----- 1. create the report directory ------------------------
+system(paste0('mkdir -p ', opt$wgcna_preprocessing_dir))
+
+
+#----- 2. generate Rmd files with Rmd templates --------------
+# a. templates without placeholder variables:
+# copy templates from tool directory to the working directory.
+# b. templates with placeholder variables:
+# substitute variables with user input values and place them in the working directory.
+
+
+#----- 01 wgcna_preprocessing.Rmd -----------------------
+readLines(opt$wgcna_preprocessing_rmd) %>%
+ (function(x) {
+ gsub('ECHO', opt$echo, x)
+ }) %>%
+ (function(x) {
+ gsub('EXPRESSION_DATA', opt$expression_data, x)
+ }) %>%
+ (function(x) {
+ gsub('OUTPUT_DIR', opt$wgcna_preprocessing_dir, x)
+ }) %>%
+ (function(x) {
+ gsub('PREPROCESSING_WORKSPACE', opt$preprocessing_workspace, x)
+ }) %>%
+ (function(x) {
+ fileConn = file('wgcna_preprocessing.Rmd')
+ writeLines(x, con=fileConn)
+ close(fileConn)
+ })
+
+
+#------ 3. render all Rmd files --------
+render('wgcna_preprocessing.Rmd', output_file = opt$wgcna_preprocessing_html)
+
+#-------4. manipulate outputs -----------------------------
+
+