generated from RobertsLab/project-template
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path02-deltadelta-boxplots copy.Rmd
193 lines (149 loc) · 6.79 KB
/
02-deltadelta-boxplots copy.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
---
title: "02-deltadelta-boxplots"
author: "Zach Bengtsson"
date: "3/22/2023"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r import data, echo=FALSE}
### import the compiled Cq data
qpcr_data <- read.csv("/Users/zachbengtsson/Desktop/Mussel Heat qPCR/Cq_Compiled_v2.csv", header = TRUE)
```
```{r delta Ct, echo=FALSE}
### calc delta Ct (GOI - NG)
qpcr_data$delta_ct_HSP24 <- qpcr_data$HSP24CqMean - 38
qpcr_data$delta_ct_HSC70 <- qpcr_data$HSC70CqMean - 38
qpcr_data$delta_ct_HSP90 <- qpcr_data$HSP90CqMean - 38
```
```{r delta delta Ct, echo=FALSE}
### Calc delta delta Ct (Subtracting Average delta Ct of control from each)
# Calculate the average delta Ct value for the control group
control_mean_delta_ct_HSP24 <- mean(qpcr_data[qpcr_data$Biological.Set.Name == "control", "delta_ct_HSP24"])
control_mean_delta_ct_HSC70 <- mean(qpcr_data[qpcr_data$Biological.Set.Name == "control", "delta_ct_HSC70"])
control_mean_delta_ct_HSP90 <- mean(qpcr_data[qpcr_data$Biological.Set.Name == "control", "delta_ct_HSP90"])
# Calculate the delta delta Ct value for each sample
qpcr_data$delta_delta_ct_HSP24 <- qpcr_data$delta_ct_HSP24 - control_mean_delta_ct_HSP24
qpcr_data$delta_delta_ct_HSC70 <- qpcr_data$delta_ct_HSC70 - control_mean_delta_ct_HSC70
qpcr_data$delta_delta_ct_HSP90 <- qpcr_data$delta_ct_HSP90 - control_mean_delta_ct_HSP90
```
```{r fold change, echo=FALSE}
### Calc the fold change (2^(-delta delta Ct))
qpcr_data$fold_change_HSP24 <- 2^(-qpcr_data$delta_delta_ct_HSP24)
qpcr_data$fold_change_HSC70 <- 2^(-qpcr_data$delta_delta_ct_HSC70)
qpcr_data$fold_change_HSP90 <- 2^(-qpcr_data$delta_delta_ct_HSP90)
```
```{r load packages, echo=FALSE}
library(ggplot2)
library(tidyr)
library(dplyr)
```
```{r data frame, echo=FALSE}
### Create data frame for boxplots
fold_change_data <- data.frame(
group = qpcr_data$Biological.Set.Name,
fold_change_HSP24 = qpcr_data$fold_change_HSP24,
fold_change_HSC70 = qpcr_data$fold_change_HSC70,
fold_change_HSP90 = qpcr_data$fold_change_HSP90
)
```
```{r reshape data, echo=FALSE}
### Reshape with tidy gather() function
fold_change_data <- fold_change_data %>%
gather(key = "gene", value = "fold_change", -group)
```
```{r boxplot, echo=FALSE}
ggplot(fold_change_data, aes(x = group, y = fold_change)) +
geom_boxplot(fill = "transparent", color = "black") +
facet_wrap(~ gene, scales = "free_y", strip.position = "top",
labeller = labeller(gene = c(fold_change_HSC70 = "HSC70", fold_change_HSP24 = "HSP24", fold_change_HSP90 = "HSP90"))) +
xlab("Group") +
ylab("Fold Change") +
ggtitle("Fold Change of Gene Expression") +
theme_bw()
ggsave("boxplot_withOutliers.png", width = 7, height = 3, dpi = 500)
```
```{r boxplot, echo=FALSE}
# Calculate upper and lower whiskers for each gene and group
whiskers <- fold_change_data %>%
group_by(group, gene) %>%
summarize(
lower_whisker = boxplot.stats(fold_change)$stats[1],
upper_whisker = boxplot.stats(fold_change)$stats[5]
)
### Whiskers for each gene and group are calculated using the boxplot.stats() function. Upper whisker is defined as the maximum value that is less than or equal to 1.5 times the interquartile range (IQR) above the third quartile, and the lower whisker is defined as the minimum value that is greater than or equal to 1.5 times the IQR below the first quartile.
# Create the plot and exclude outliers
ggplot(fold_change_data, aes(x = reorder(group, fold_change, median), y = fold_change)) +
geom_boxplot(aes(color = group), fill = "transparent") +
facet_wrap(~ gene, scales = "free_y", strip.position = "top",
labeller = labeller(gene = c(fold_change_HSC70 = "HSC70", fold_change_HSP24 = "HSP24", fold_change_HSP90 = "HSP90"))) +
xlab("Group") +
ylab("Fold Change") +
ggtitle("Fold Change of Gene Expression") +
theme_bw() +
scale_color_manual(values = c("variable" = "#548235","chronic" = "#C55A11", "control" = "#2F5597")) +
scale_y_continuous(limits = c(min(whiskers$lower_whisker), max(whiskers$upper_whisker)))
ggsave("boxplot_removeOutliers.png")
```
```{r, echo=FALSE}
ggplot(fold_change_data, aes(x = reorder(group, fold_change, median), y = fold_change)) +
geom_boxplot(aes(fill = group), color = "black") +
facet_wrap(~ gene, scales = "free_y", strip.position = "top",
labeller = labeller(gene = c(fold_change_HSC70 = "HSC70", fold_change_HSP24 = "HSP24", fold_change_HSP90 = "HSP90"))) +
xlab("Group") +
ylab("Fold Change") +
ggtitle("Fold Change of Gene Expression") +
theme_bw() +
scale_fill_manual(values = c("chronic" = "blue", "variable" = "green", "control" = "red")) +
guides(fill=guide_legend(title="Group")) +
scale_y_continuous(limits = c(min(fold_change_data$fold_change), max(fold_change_data$fold_change)))
```
```{r split data, echo=FALSE}
# create separate data frames for each gene
HSC70_fold <- subset(fold_change_data, select = c("group", "fold_change"), subset = gene == "fold_change_HSC70")
HSP90_fold <- subset(fold_change_data, select = c("group", "fold_change"), subset = gene == "fold_change_HSP90")
HSP24_fold <- subset(fold_change_data, select = c("group", "fold_change"), subset = gene == "fold_change_HSP24")
```
```{r ANOVA for each gene, echo=FALSE}
# run ANOVA for each gene data frame
HSC70_aov <- aov(fold_change ~ group, data = HSC70_fold)
HSP90_aov <- aov(fold_change ~ group, data = HSP90_fold)
HSP24_aov <- aov(fold_change ~ group, data = HSP24_fold)
# display ANOVA summary
summary(HSC70_aov)
summary(HSP90_aov)
summary(HSP24_aov)
# perform Tukey's HSD test for each gene data frame
HSC70_tukey <- TukeyHSD(HSC70_aov)
HSP90_tukey <- TukeyHSD(HSP90_aov)
HSP24_tukey <- TukeyHSD(HSP24_aov)
# display pairwise comparison results
print(HSC70_tukey)
print(HSP90_tukey)
print(HSP24_tukey)
```
```{r ANOVA no outliers, echo=FALSE}
# remove outliers from HSC70_fold
HSC70_fold_no_outliers <- HSC70_fold[!HSC70_fold$fold_change %in% boxplot.stats(HSC70_fold$fold_change)$out,]
# remove outliers from HSP90_fold
HSP90_fold_no_outliers <- HSP90_fold[!HSP90_fold$fold_change %in% boxplot.stats(HSP90_fold$fold_change)$out,]
# remove outliers from HSP24_fold
HSP24_fold_no_outliers <- HSP24_fold[!HSP24_fold$fold_change %in% boxplot.stats(HSP24_fold$fold_change)$out,]
# run ANOVA for each gene data frame
HSC70_no_aov <- aov(fold_change ~ group, data = HSC70_fold_no_outliers)
HSP90_no_aov <- aov(fold_change ~ group, data = HSP90_fold_no_outliers)
HSP24_no_aov <- aov(fold_change ~ group, data = HSP24_fold_no_outliers)
# display ANOVA summary
summary(HSC70_no_aov)
summary(HSP90_no_aov)
summary(HSP24_no_aov)
# perform Tukey's HSD test for each gene data frame
HSC70_no_tukey <- TukeyHSD(HSC70_no_aov)
HSP90_no_tukey <- TukeyHSD(HSP90_no_aov)
HSP24_no_tukey <- TukeyHSD(HSP24_no_aov)
# display pairwise comparison results
print(HSC70_no_tukey)
print(HSP90_no_tukey)
print(HSP24_no_tukey)
```