Mercurial > repos > mingchen0919 > rmarkdown_deseq2

comparison DESeq_results.Rmd @ 6:2f8ddef8d545 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update deseq2
author mingchen0919
date Tue, 07 Nov 2017 13:50:32 -0500
parents 7231d7e8d3ed
children
comparison
equal deleted inserted replaced
5:fd3514267506 6:2f8ddef8d545
8 highlight: tango 8 highlight: tango
9 --- 9 ---
10 10
11 ```{r setup, include=FALSE, warning=FALSE, message=FALSE} 11 ```{r setup, include=FALSE, warning=FALSE, message=FALSE}
12 knitr::opts_chunk$set( 12 knitr::opts_chunk$set(
13 echo = ECHO 13 echo = ECHO,
14 error = TRUE
14 ) 15 )
15
16 library(DESeq2)
17 library(pheatmap)
18 library(genefilter)
19 ``` 16 ```
20 17
21 # Import workspace
22 18
23 ```{r eval=TRUE} 19 ```{r eval=TRUE}
20 # Import workspace
24 fcp = file.copy("DESEQ_WORKSPACE", "deseq.RData") 21 fcp = file.copy("DESEQ_WORKSPACE", "deseq.RData")
25 load("deseq.RData") 22 load("deseq.RData")
26 ``` 23 ```
27 24
28 # Results {.tabset} 25 # Results {.tabset}
29 26
30 ## Result table 27 ## Result table
31 28
32 ```{r} 29 ```{r}
33 group = colnames(sample_table)[CONTRAST_GROUP] 30 cat('--- View the top 100 rows of the result table ---')
34 res <- results(dds, contrast = c(group, 'TREATMENT_LEVEL', 'CONDITION_LEVEL')) 31 res <- results(dds, contrast = c('CONTRAST_FACTOR', 'TREATMENT_LEVEL', 'CONDITION_LEVEL'))
35 datatable(as.data.frame(res), style="bootstrap", filter = 'top', 32 write.csv(as.data.frame(res), file = 'deseq_results.csv')
33 res_df = as.data.frame(res)[1:100, ]
34 datatable(res_df, style="bootstrap", filter = 'top',
36 class="table-condensed", options = list(dom = 'tp', scrollX = TRUE)) 35 class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
37 ``` 36 ```
38 37
39 ## Result summary 38 ## Result summary
40 39
43 ``` 42 ```
44 43
45 44
46 # MA-plot {.tabset} 45 # MA-plot {.tabset}
47 46
48 ## Shrinked with `lfcShrink()` function
49 47
50 ```{r eval=FALSE}
51 shrink_res = DESeq2::lfcShrink(dds, contrast = c(group, 'TREATMENT_LEVEL', 'CONDITION_LEVEL'), res=res)
52 plotMA(shrink_res)
53 ```
54
55 ## Shrinked with Bayesian procedure
56 48
57 ```{r} 49 ```{r}
50 cat('--- Shrinked with Bayesian procedure ---')
58 plotMA(res) 51 plotMA(res)
59 ``` 52 ```
60 53
61 54
62 # Histogram of p values 55 # Histogram of p values
66 col = "grey50", border = "white", main = "", 59 col = "grey50", border = "white", main = "",
67 xlab = "Mean normalized count larger than 1") 60 xlab = "Mean normalized count larger than 1")
68 ``` 61 ```
69 62
70 63
71 # Gene clustering 64 # Visualization {.tabset}
65 ## Gene clustering
72 66
73 ```{r} 67 ```{r}
74 group_index = as.numeric(strsplit("CLUSTERING_GROUPS", ',')[[1]]) 68 clustering_groups = strsplit("CLUSTERING_FACTORS", ',')[[1]]
75 clustering_groups = colnames(sample_table)[group_index]
76 69
77 topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20) 70 topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
78 mat <- assay(rld)[ topVarGenes, ] 71 mat <- assay(rld)[ topVarGenes, ]
79 mat <- mat - rowMeans(mat) 72 mat <- mat - rowMeans(mat)
80 annotation_col <- as.data.frame(colData(rld)[, clustering_groups]) 73 annotation_col <- as.data.frame(colData(rld)[, clustering_groups])
81 colnames(annotation_col) = clustering_groups 74 colnames(annotation_col) = clustering_groups
82 rownames(annotation_col) = colnames(mat) 75 rownames(annotation_col) = colnames(mat)
83 pheatmap(mat, annotation_col = annotation_col) 76 pheatmap(mat, annotation_col = annotation_col)
84 ``` 77 ```
85 78
79 ## Sample-to-sample distance
80
81 ```{r}
82 sampleDistMatrix <- as.matrix( sampleDists )
83 colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
84 pheatmap(sampleDistMatrix,
85 clustering_distance_cols = sampleDists,
86 col = colors)
87 ```
88
89 ## PCA plot
90
91 ```{r}
92 plotPCA(rld, intgroup = clustering_groups)
93 ```
94
95 ## MDS plot {.tabset}
96
97 ### Data table
98 ```{r}
99 mds <- as.data.frame(colData(rld)) %>%
100 cbind(cmdscale(sampleDistMatrix))
101 knitr::kable(mds)
102 ```
103
104 ### Plot
105 ```{r}
106 ggplot(mds, aes(x = `1`, y = `2`, col = time)) +
107 geom_point(size = 3) + coord_fixed()
108 ```
109