%\VignetteEngine{knitr::knitr}
\documentclass[a4paper]{article}

\title{Usage of MODA}
\author{Dong Li \\ dxl466@cs.bham.ac.uk\\
School of Computer Science, The University of Birmingham, UK
}
\date{Date modified: \Sexpr{Sys.Date()}}
%\VignetteIndexEntry{Usage of MODA}
<<style-Sweave, eval=TRUE, echo=FALSE, results='asis'>>=
BiocStyle::latex()
@
\begin{document}

<<include=FALSE>>=
library(knitr)
opts_chunk$set(
concordance=TRUE
)
@

\maketitle

In this example we embed parts of the examples from the
\texttt{MODA} help page into a single document.

\section{Module detection}
First of all we conduct the experiment on the synthetic dataset which contains two expression profiles $datExpr1$ and $datExpr2$ with 500 genes, and each has 20 and 25 samples. Details of data generation can be found in supplementary file of MODA paper \cite{MODA}. Basic module detection functions are provided by WGCNA \cite{WGCNA}.

<<>>=
library(MODA)
data(synthetic)
ResultFolder = 'ForSynthetic' # where middle files are stored
CuttingCriterion = 'Density' # could be Density or Modularity
indicator1 = 'X'     # indicator for data profile 1
indicator2 = 'Y'      # indicator for data profile 2
specificTheta = 0.1 #threshold to define condition specific modules
conservedTheta = 0.1#threshold to define conserved modules
##modules detection for network 1
intModules1 <- WeightedModulePartitionDensity(datExpr1,ResultFolder,
                                indicator1,CuttingCriterion)
##modules detection for network 2
intModules2 <- WeightedModulePartitionDensity(datExpr2,ResultFolder,
                                indicator2,CuttingCriterion)
@
which shows how to detect modules using hierecal clustering with the optimal cutting height of dendrogram. The heatmap of correlation matrix of gene expression profile 1 may looks like Figure \ref{fig:heat1}. Another package \cite{link} has the similar function.
\begin{figure}
\centering
<<echo=FALSE,out.width='0.6\\textwidth'>>=
heatmap(cor(as.matrix(datExpr1)))
@
\caption{Correlation matrix of gene expression profile 1.}
\label{fig:heat1}
\end{figure}

The selection of optimal cutting height for each expression profile would be stored under directory $ResultFolder$. Take $datExpr1$ in the synthetic data for example, a file named $Partitions\_X.pdf$ may looks like Figure \ref{fig:cut1}.
\begin{figure}
\centering
\includegraphics[width=10cm]{ForSynthetic/Partitions_X.pdf}
\caption{Maximal partition density based hierarchical clustering for gene expression profile 1.}
\label{fig:cut1}
\end{figure}

At the same time, each module for each expression profile would be stored as plain text file, with the name indicator from $indicator1$ and $indicator2$. Each secondary directory under $ResultFolder$ has the same name of condition name, e.g $indicator2$, used to store differential analysis results.

\section{Network comparison}
After the module detection for background network and all condition-specific networks, we can compare them using following function
<<>>=
CompareAllNets(ResultFolder,intModules1,indicator1,intModules2,
               indicator2,specificTheta,conservedTheta)
@

The condition specific networks can be specified by two vectors if there are more. There are three files under the secondary directory named by condition name: two text files of them are condition specific and conserved modules id from background network, and one pdf for showing how to determine these modules by two parameters $specificTheta$ and $conservedTheta$ based on a Jaccard index matrix. Theoretical details can be found in supplementary file of MODA paper. The figure may looks like Figure \ref{fig:compare}.
\begin{figure}
\centering
\includegraphics[width=10cm]{ForSynthetic/Y/module_overlap_removeY.pdf}
\caption{Overlap degree of modules from two networks.}
\label{fig:compare}
\end{figure}

\section{Biological explanation}
Finally we can do gene annotation enrichment analysis with intergative tools like DAVID\footnote{https://david.ncifcrf.gov} or Enrichr\footnote{http://amp.pharm.mssm.edu/Enrichr}, to see whether a module gene list can be explained by existing biological process, pathways or even diseases.
\section{Session info}
<<sessionInfo, results='asis', echo=FALSE>>=
toLatex(sessionInfo())
@
 \begin{thebibliography}{1}

  \bibitem{MODA} Dong Li et al. {\em MODA: MOdule Differential Analysis for weighted gene co-expression network} bioRxiv 053496 (2016).
  
  \bibitem{WGCNA} Langfelder, Peter, and Steve Horvath. {\em WGCNA: an R package for weighted correlation network analysis.} BMC bioinformatics 9.1 (2008): 1.
  
  \bibitem{link} Kalinka, Alex T., and Pavel Tomancak. {\em linkcomm: an R package for the generation, visualization, and analysis of link communities in networks of arbitrary size and type.} Bioinformatics 27.14 (2011).

  \end{thebibliography}
\end{document}