Mercurial > repos > mingchen0919 > rmarkdown_wgcna

comparison wgcna_eigengene_visualization.Rmd @ 0:4275479ada3a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/statonlab/docker-GRReport/tree/master/my_tools/rmarkdown_wgcna commit d91f269e8bc09a488ed2e005122bbb4a521f44a0-dirty
author mingchen0919
date Tue, 08 Aug 2017 12:35:50 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:4275479ada3a
1 ---
2 title: 'WGCNA: eigengene visualization'
3 output:
4 html_document:
5 number_sections: true
6 toc: true
7 theme: cosmo
8 highlight: tango
9 ---
10
11 ```{r setup, include=FALSE, warning=FALSE, message=FALSE}
12 knitr::opts_chunk$set(
13 echo = ECHO
14 )
15 ```
16
17 # Import workspace
18
19 This step imports workspace from the **WGCNA: construct network** step.
20
21 ```{r}
22 fcp = file.copy("CONSTRUCT_NETWORK_WORKSPACE", "deseq.RData")
23 load("deseq.RData")
24 ```
25
26
27 # Gene modules {.tabset}
28
29 ```{r}
30 if(!is.na(SOFT_THRESHOLD_POWER)) soft_threshold_power = SOFT_THRESHOLD_POWER
31 ```
32
33 ## Identify gene modules
34
35 The gene network is constructed based on **soft threshold power = `r soft_threshold_power`**
36
37 ```{r}
38 gene_network = blockwiseModules(expression_data, power = soft_threshold_power,
39 TOMType = "unsigned", minModuleSize = 30,
40 reassignThreshold = 0, mergeCutHeight = 0.25,
41 numericLabels = TRUE, pamRespectsDendro = FALSE,
42 verbose = 3)
43 ```
44
45
46 ```{r}
47 modules = table(gene_network$colors)
48 n_modules = length(modules) - 1
49 module_size_upper = modules[2]
50 module_size_lower = modules[length(modules)]
51
52 module_table = data.frame(model_label = c(0, 1:n_modules),
53 gene_size = as.vector(modules))
54 datatable(t(module_table))
55 ```
56
57 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.
58
59
60 ## Dendrogram and module plot
61
62 ```{r}
63 # Convert labels to colors for plotting
64 module_colors = labels2colors(gene_network$colors)
65 # Plot the dendrogram and the module colors underneath
66 plotDendroAndColors(gene_network$dendrograms[[1]], module_colors[gene_network$blockGenes[[1]]],
67 "Module colors",
68 dendroLabels = FALSE, hang = 0.03,
69 addGuide = TRUE, guideHang = 0.05)
70 ```
71
72
73 # Gene module correlation
74
75 We can calculate eigengenes and use them as representative profiles to quantify similarity of found gene modules.
76
77 ```{r}
78 n_genes = ncol(expression_data)
79 n_samples = nrow(expression_data)
80 ```
81
82 ```{r}
83 diss_tom = 1-TOMsimilarityFromExpr(expression_data, power = soft_threshold_power)
84 set.seed(123)
85 select_genes = sample(n_genes, size = PLOT_GENES)
86 select_diss_tom = diss_tom[select_genes, select_genes]
87
88 # calculate gene tree on selected genes
89 select_gene_tree = hclust(as.dist(select_diss_tom), method = 'average')
90 select_module_colors = module_colors[select_genes]
91
92 # transform diss_tom with a power to make moderately strong connections more visiable in the heatmap.
93 plot_diss_tom = select_diss_tom^7
94 # set diagonal to NA for a nicer plot
95 diag(plot_diss_tom) = NA
96 ```
97
98
99 ```{r fig.align='center'}
100 TOMplot(plot_diss_tom, select_gene_tree, select_module_colors, main = "Network heatmap")
101 ```
102
103
104 # Eigengene visualization {.tabset}
105
106 ## Eigengene dendrogram
107
108 ```{r fig.align='center'}
109 module_eigengenes = moduleEigengenes(expression_data, module_colors)$eigengenes
110 plotEigengeneNetworks(module_eigengenes, "Eigengene dendrogram",
111 plotHeatmaps = FALSE)
112 ```
113
114 ## Eigengene adjacency heatmap
115
116 ```{r fig.align='center'}
117 plotEigengeneNetworks(module_eigengenes, "Eigengene adjacency heatmap",
118 marHeatmap = c(2, 3, 2, 2),
119 plotDendrograms = FALSE, xLabelsAngle = 90)
120 ```
121