Establish design & contrast matrices
## Establish design matrix
design = model.matrix(~ 0 +
group +
age +
sex +
pool)
colnames(design) <- gsub("group", "", colnames(design))
contr.matrix = makeContrasts(
CtlvI1 = Ctl - I1,
CtlvI2 = Ctl - I2,
CtlvI5 = Ctl - I5,
CtlvI15 = Ctl - I15,
CtlvI29 = Ctl - I29,
I1vI2 = I1 - I2,
I1vI5 = I1 - I5,
I1vI15 = I1 - I15,
I1vI29 = I1 - I29,
I2vI5 = I2 - I5,
I2vI15 = I2 - I15,
I2vI29 = I2 - I29,
I5vI15 = I5 - I15,
I5vI29 = I5 - I29,
I15vI29 = I15 - I29,
levels = colnames(design))
## Mean-variance trend and sample weight plots
v <- voomWithQualityWeights(dge, design, plot = TRUE)

## Fitting the model
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
tfit <- treat(vfit, lfc = 0.5)
plotSA(tfit)

Volcano plot
tmp1 <- topTreat(tfit, coef = 1, n = Inf)
results1 <- mutate(tmp1, sig=ifelse(tmp1$adj.P.Val<0.1, "Sig", "Not Sig"))
p1 <- ggplot(results1, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Control vs. Infected, Day 1") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp2 <- topTreat(tfit, coef = 2, n = Inf)
results2 <- mutate(tmp2, sig=ifelse(tmp2$adj.P.Val<0.1, "Sig", "Not Sig"))
p2 <- ggplot(results2, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Control vs. Infected, Day 2") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp3 <- topTreat(tfit, coef = 3, n = Inf)
results3 <- mutate(tmp3, sig=ifelse(tmp3$adj.P.Val<0.1, "Sig", "Not Sig"))
p3 <- ggplot(results3, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Control vs. Infected, Day 5") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp4 <- topTreat(tfit, coef = 4, n = Inf)
results4 <- mutate(tmp4, sig=ifelse(tmp4$adj.P.Val<0.1, "Sig", "Not Sig"))
p4 <- ggplot(results4, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Control vs. Infected, Day 15") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp5 <- topTreat(tfit, coef = 5, n = Inf)
results5 <- mutate(tmp5, sig=ifelse(tmp5$adj.P.Val<0.1, "Sig", "Not Sig"))
p5 <- ggplot(results5, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Control vs. Infected, Day 29") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp6 <- topTreat(tfit, coef = 6, n = Inf)
results6 <- mutate(tmp6, sig=ifelse(tmp6$adj.P.Val<0.1, "Sig", "Not Sig"))
p6 <- ggplot(results6, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 1 vs. Infected, Day 2") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp7 <- topTreat(tfit, coef = 7, n = Inf)
results7 <- mutate(tmp7, sig=ifelse(tmp7$adj.P.Val<0.1, "Sig", "Not Sig"))
p7 <- ggplot(results7, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 1 vs. Infected, Day 5") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp8 <- topTreat(tfit, coef = 8, n = Inf)
results8 <- mutate(tmp8, sig=ifelse(tmp8$adj.P.Val<0.1, "Sig", "Not Sig"))
p8 <- ggplot(results8, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 1 vs. Infected, Day 15") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp9 <- topTreat(tfit, coef = 9, n = Inf)
results9 <- mutate(tmp9, sig=ifelse(tmp9$adj.P.Val<0.1, "Sig", "Not Sig"))
p9 <- ggplot(results9, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 1 vs. Infected, Day 29") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp10 <- topTreat(tfit, coef = 10, n = Inf)
results10 <- mutate(tmp10, sig=ifelse(tmp10$adj.P.Val<0.1, "Sig", "Not Sig"))
p10 <- ggplot(results10, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 2 vs. Infected, Day 5") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp11 <- topTreat(tfit, coef = 11, n = Inf)
results11 <- mutate(tmp11, sig=ifelse(tmp11$adj.P.Val<0.1, "Sig", "Not Sig"))
p11 <- ggplot(results11, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 2 vs. Infected, Day 15") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp12 <- topTreat(tfit, coef = 12, n = Inf)
results12 <- mutate(tmp12, sig=ifelse(tmp12$adj.P.Val<0.1, "Sig", "Not Sig"))
p12 <- ggplot(results12, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 2 vs. Infected, Day 29") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp13 <- topTreat(tfit, coef = 13, n = Inf)
results13 <- mutate(tmp13, sig=ifelse(tmp13$adj.P.Val<0.1, "Sig", "Not Sig"))
p13 <- ggplot(results13, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 5 vs. Infected, Day 15") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp14 <- topTreat(tfit, coef = 14, n = Inf)
results14 <- mutate(tmp14, sig=ifelse(tmp14$adj.P.Val<0.1, "Sig", "Not Sig"))
p14 <- ggplot(results14, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 5 vs. Infected, Day 29") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
tmp15 <- topTreat(tfit, coef = 15, n = Inf)
results15 <- mutate(tmp15, sig=ifelse(tmp15$adj.P.Val<0.1, "Sig", "Not Sig"))
p15 <- ggplot(results15, aes(logFC, -log10(P.Value))) +
geom_point(aes(col = sig)) +
scale_color_manual(values=c("black", "red")) +
ggtitle("Infected, Day 15 vs. Infected, Day 29") +
ylab("-Log10(p-value)") +
xlab("Log2(Fold change)") +
theme_bw() +
theme(legend.position = "none")
grid.arrange(p1, p2, p3, p4, p5, p6, p7,
p8, p9, p10, p11, p12, p13, p14, p15, nrow = 8)

Boxplots for differentially expressed genes
lcpm.DE <- lcpm %>%
as_tibble(rownames = NA) %>%
rownames_to_column("identifier") %>%
filter(identifier %in% dt.tib$gene) %>%
pivot_longer(cols = contains("_"),
names_to = "sample",
values_to = "lcpm") %>%
separate(sample, into = c("tissue", "bird")) %>%
mutate(bird = as.factor(bird)) %>%
left_join(covars, by = "bird") %>%
select(identifier, bird, lcpm, group) %>%
mutate(group = recode(group, C29 = "Control"),
group = recode(group, C1 = "Control"),
tmpID = ifelse(group == "Control", 'Control', 'Infected'))
aprioriPlotting <- function(target, ...) {
annot.target <- annot %>%
select(
ensembl_gene_id,
hgnc_symbol) %>%
filter(ensembl_gene_id == target) %>%
mutate(hgnc_symbol = replace_na(hgnc_symbol, "."))
plot <- lcpm.DE %>%
filter(identifier == target) %>%
mutate(group = fct_relevel(group, "Control", "I1", "I2", "I5", "I15", "I29")) %>%
ggplot(aes(x = group, y = lcpm)) +
facet_grid(. ~ tmpID, scales = "free", space = "free") +
ylab("Log2(Counts per million)") +
xlab("Group") +
geom_point(position = position_dodge(width=0.75), show.legend = FALSE) +
geom_boxplot(alpha = 0.5) +
geom_label_repel(aes(label = bird, fill = NULL), show.legend=FALSE) +
theme_classic() +
labs(title= paste0(target, " - ", annot.target[1,2])) +
theme(legend.title = element_blank())
print(plot)
}
for(target in sort(unique(lcpm.DE$identifier))) {
aprioriPlotting(target)
}


