-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathanalysis8hpaAllR.Rmd
executable file
·244 lines (165 loc) · 10.6 KB
/
analysis8hpaAllR.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
---
title: "8hpa Control - All Replicates"
author: "Kirsten Gotting"
output:
html_document:
code_download: true
toc: true
toc_depth: 2
toc_float:
collapsed: false
---
```{r initializeData8hpaControlAllReplicates, echo = FALSE, results = "hide", message = FALSE, error = FALSE, warning = FALSE, cache = FALSE}
## user input to 'DESeq2All.Rdata'
load("./output/DESeq2AllAllR.Rdata")
setwd("./")
## First list the R libraries I want to use
libraries <- c("DESeq2", "gplots", "RColorBrewer", 'knitr', 'grid', 'gridExtra', "tidyverse", 'kiRsten')
## Now read in the libraries
lapply(libraries, function(x){
library(x, character.only = TRUE, quietly = TRUE)
})
## Initialize knitr options
opts_chunk$set(echo=FALSE, message=FALSE, results="hide", fig.keep="all", warning=FALSE, error=FALSE, fig.path="figures/", cache = TRUE)
opts_knit$set(root.dir = "./")
```
# Overview
Differential expression analysis with 8hpa control and all replicates.
#Spearman Correlation
Spearman's rank correlation coefficient is a nonparametric measure of rank correlation, recommended for data that does not necessarily come from a bivariate normal distribution (Wikipedia). Below, we plot the Spearman correlation of the RPM (reads per million) for all genes that have an average RPM ≥ 1.
```{r spearman_correlation8hpaControlAllReplicates, fig.cap = ""}
## Make a spearman plot
filteredRPMs <- normCounts %>% group_by(ID) %>% mutate(rpmFlag = ifelse((mean(rpm) >= 1), yes = 'yes', no = 'no')) %>% ungroup() %>%
filter(rpmFlag == 'yes') %>% select(ID, Sample, rpm) %>% spread(Sample, rpm)
rpmfilt <- data.frame(select(filteredRPMs, -ID), row.names = filteredRPMs$ID, check.names = FALSE)
spearman <- data.frame(cor(rpmfilt, method = 'spearman'), check.names = FALSE)
spearman$ID <- rownames(spearman)
spearmanTidy <- spearman %>% gather(Sample, correlation, -ID)
## and plot it!
ggplot(data = spearmanTidy, aes(x = ID, y = Sample, fill = correlation)) + geom_tile() +
scale_fill_gradient2(low = "yellow", mid = 'white', high = "purple4", midpoint = 0.8, limit = c(0,1)) +
ggtitle('Spearman Correlation of Samples') + xlab('') + ylab('') +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
```
# Pairwise PCA
Principal component analysis (PCA) is a statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components (Wikipedia). Below is a PCA of the 8hpf control versus the other developmental timepoints.
```{r PCA8hpaControlAllReplicates, fig.fullwidth = TRUE, fig.width = 10, fig.height = 10}
pcas_data <- lapply(contrast_names, function(x){
y <- plotPCA(rld_all[[x]], returnData = TRUE, ntop = 10000)
y$ID <- x
y
})
pca_data <- bind_rows(pcas_data)
ggplot(data = pca_data, aes(x = PC1, y = PC2, label = name, colour = factor(condition))) +
geom_point(size = 3) + theme_bw() + geom_text(hjust = 0, vjust = -0.4) +
xlab("PC 1") + ylab("PC 2") +
scale_x_continuous(expand = c(0.98, 0)) +
scale_y_continuous(expand = c(0.2, 0)) +
theme(legend.position="none") +
facet_wrap(~ ID, scales = "free")
```
# MA and Volcano plots {.tabset}
MA plot: The log2 fold change (M) plotted against the log2 average (A) of the normalized read count for each gene. Significantly differentially expressed genes are plotted in red.
Volcano plot: The log2 fold change (M) plotted against the -log10 (eg. 1e^-10 = 10) of the adjusted p-value.
```{r mavol8hpaControlAllReplicates}
x <- contrast_names[[1]]
plots <- lapply(contrast_names, function(x){
df <- all_results_tidyDF %>% filter(contrastID == x) %>% spread(dea_ID, dea_Value) %>% mutate(logbaseMean = log(baseMean))
pval <- as.numeric(pval)
if(is.null(df)){
message(paste0("Pvalue too stringent for ", contrast, ". MA/Volcano plot ommited.\n"))
} else {
pval <- as.numeric(pval)
lower_label <- paste0("Genes with Counts: ", nrow(df), "\nUp: ",#create the label of the legend
nrow(subset(df, log2FoldChange > 0 & padj < pval)),
", Down: ", nrow(subset(df, log2FoldChange < 0 & padj < pval)), "\n")
df <- df %>% mutate(significance_group = ifelse(padj <= pval, yes = 'yes', no = 'no'))
legend_labels <- c(paste0("padj > ", pval,": ", nrow(filter(df, significance_group == 'no'))),
paste0("padj < ", pval,": ", nrow(filter(df, significance_group == 'yes'))))
ma.plot <- ggplot(df, aes(x = logbaseMean, y = log2FoldChange, group = significance_group, colour = significance_group)) +
geom_point() + theme_bw() +
ggtitle('MA plot') + scale_x_continuous(expand = c(0, 0)) +
xlab("A = Mean of log2(Abundance)") + ylab("M = log2(Fold Change)") +
scale_color_manual(values = c("gray24", "red"), name = lower_label, breaks = c("gray24", "red"), labels = legend_labels) +
geom_hline(aes(yintercept = 0), colour = "black", linetype = "dashed") +
geom_hline(aes(yintercept = 1), colour = "gray66", linetype = "dashed") +
geom_hline(aes(yintercept = -1), colour = "gray66", linetype="dashed") +
theme(legend.position = 'none', text = element_text(size=20))
volcano.plot <- ggplot(df, aes(x = log2FoldChange, y = -log(padj, base = c(10)), group = significance_group, colour = significance_group)) +
geom_point() +
ggtitle('Volcano plot') +
xlab("M = log2(Fold Change)") + ylab("-log10(Adjusted P-Value)") +
scale_color_manual(values = c("gray24", "red"), labels = legend_labels, name=lower_label) +
theme_bw() +
theme(legend.title = element_text(size = 16, face = "bold"), legend.text = element_text(size = 20), text = element_text(size=20))
grob <- marrangeGrob(list(ma.plot,
volcano.plot),
top = textGrob(x, gp = gpar(fonsize = 40, fontface = 'bold', cex = 2), hjust = 1.8),
nrow = 1, ncol = 2,
widths=c(0.38, 0.62))
invisible(grob)
}
})
names(plots) <- contrast_names
```
```{r plotmavol8hpaControlAllReplicates, fig.keep = "all", fig.height = 10, fig.width = 20, echo=FALSE, results='asis'}
for(i in contrast_names){
cat(paste0("\n\n\n##", i, "\n"))
print(plots[[i]])
}
```
#Heatmaps
The first heatmap was made by hierarchically clustering the euclidean distances of the log2 fold change for genes with p-values less than 1e^-5 at any time point.
The second heatmap uses the same hierarchical method to cluster timepoints as well as genes.
```{r pairwiseHeatmaps8hpaControlAllReplicates, results = "hide", fig.cap = "", fig.keep = "high"}
lfc_table <- all_results_tidyDF %>%
mutate(significant = ifelse(dea_ID == 'padj' & dea_Value <= pval, yes = TRUE, no = FALSE)) %>% # label all significant genes
group_by(contrastID, Gene) %>% # label the # of times genes that are significant
ungroup() %>%
group_by(Gene) %>%
filter(any(significant), dea_ID == 'log2FoldChange') %>%
select(-significant) %>% unite(colname, contrastID, dea_ID) %>% # make the table wide formatted
spread(colname, dea_Value)
## create a table of significant genes with the cluster colors assigned
sign.table <- all_results_tidyDF %>%
spread(key = dea_ID, value = dea_Value) %>%
filter(Gene %in% lfc_table$Gene) %>%
group_by(contrastID, Gene) %>%
mutate(sortby = -log(padj, base = c(10))*sign(log2FoldChange)) %>%
ungroup() %>% gather(dea_ID, dea_Value, -Gene, -contrastID) %>%
unite(idAll, contrastID, dea_ID) %>% spread(key = idAll, value = dea_Value) %>%
dplyr::rename('TranscriptID' = Gene)
heatmap.input.table <- data.frame(lfc_table)
rownames(heatmap.input.table) <- lfc_table$Gene
heatmap.input.table$Gene <- NULL
heatmap.input.table <- na.omit(data.matrix(heatmap.input.table[which(rownames(heatmap.input.table) %in% sign.table$TranscriptID), ]))
colnames(heatmap.input.table) <- contrast_names
prettycolors <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
col.breaks <- seq(-3, 3, length.out = ncol(heatmap.input.table))
## Make two heatmaps, allow custimizable graphing
hr <- hclust(dist(heatmap.input.table), method="complete") # Creates row dendrogram
hc <- hclust(dist(t(heatmap.input.table)), method="complete") # Creates column dendrogram
myclr <- cutree(hr, k = 7) # set up where the groups will be, change 'k' for more or less groups
mycolr <- c("Green", "Red", "Blue", "Orange", "Purple", "Pink", "Yellow") # must be the same length as 'k'
mycolr <- mycolr[as.vector(myclr)] # pull out the colors assigned to each gene
mycolr.df <- data.frame(Cluster_Color = mycolr, TranscriptID = names(myclr)) # create a data frame of the cluster color each gene belongs to
if(ncol(heatmap.input.table) > 1){
## Make a heatmap without a column dendrogram
col.breaks <- seq(-3, 3, length.out = ncol(heatmap.input.table))
heatmap <- heatmap.2(heatmap.input.table, margin = c(15,5), col = prettycolors, dendrogram = 'row', trace = c("none"), Rowv=as.dendrogram(hr), Colv = NA, labRow = FALSE, key = TRUE, breaks = col.breaks, key.ylab = "Gene Count", key.xlab = "log2FoldChange", cexCol = 1.5, RowSideColors = mycolr)
title(main = paste0("Clustering of ", nrow(heatmap.input.table)," Genes"), line = 0)
par(cex.main=0.9)
}
if(ncol(heatmap.input.table) > 2) {
k <- ifelse(ncol(heatmap.input.table < 5), 3, 5)
myclc <- cutree(hc, k = k) # create groupings for the column dendrogram
mycolc <- c("Green", "Red", "Blue", "Orange", "Purple")
mycolc <- mycolc[as.vector(myclc)]
## Make a heatmap of the logfold changes with column dendrogram.
heatmap <- heatmap.2(heatmap.input.table, margin = c(15,5), dendrogram = 'both', Rowv = as.dendrogram(hr), Colv = as.dendrogram(hc), col = prettycolors, trace = c("none"), labRow = FALSE, key = TRUE, breaks = col.breaks, main = paste0("Clustering of Timepoints on ", nrow(heatmap.input.table)," Genes"), key.ylab = "Gene Count", key.xlab = "log2FoldChange", cexCol = 1.5, RowSideColors = mycolr, ColSideColors = mycolc)
}
```
```{r session8hpaControlAllReplicates, echo=FALSE}
writeLines(capture.output(sessionInfo()), "sessionInfo/8hpaAllSessionInfo.txt")
```
#[R-session information](sessionInfo/8hpaAllSessionInfo.txt)