Comparing standard runs

note

Thius script was updated on Jan 24,2020, to exclude the lowest dilution series. The original data was saved as “two.plate.standards.backup.csv”

The goals of this script are to 1. create a simple xy-plot of two standard qpcr runs, consisting each of ~150 primers/probes run at eight different concentrations in triplicate.
* The points need to be color-coded by template DNA concentration.
* linear correlation line and stats * 95% confidence interval * 95% prediction interval * 95% confidence scatterplot-based ellipse

create boxplots for each concentration for each plate

Figures with various settings are displayed below… all figures have a linear regression line with 99% confidence interval. NOTE!!! : The confidence interval is tiny, but you can zoom in to be sure that it is there

raw.data <- read.csv("two.plate.standards.csv", header = T)
raw.data$copies.33nl.well <- as.numeric(raw.data$copies.33nl.well)
raw.data$copies.33nl.well <- formatC(raw.data$copies.33nl.well, format = "e", digits = 2)

Figure 1

Full-data confidence ellipse

library(ggplot2)
library(scales)

my.color = "#6494e0"

df = raw.data
df = df[complete.cases(df),]
#df = df[df$Target_Name=="BTUC_GTVS",]
#options(scipen=999)
theme_set(theme_bw())
#lreg.gglot <- abline(lm(raw.data$mix_1_Ct ~ raw.data$mix_2_Ct))
plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) +xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), level = 0.95, size = 0.5, color = "#4c4c4c",linetype = 1) +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) +
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 1") + guides(color = guide_legend(override.aes = list(size=5)))+
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
#confidenceEllipse(lm(raw.data$mix_1_Ct ~ raw.data$mix_2_Ct), add = T, levels = 0.95)
#geom_line(lreg.ggplot)
ggsave(file="2ctA.pdf", device="pdf", plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

Figure 2

By-concentration confidence ellipse

library(ggplot2)
library(scales)

my.color = "#6494e0"

df = raw.data
df = df[complete.cases(df),]
#df = df[df$Target_Name=="BTUC_GTVS",]
#options(scipen=999)
theme_set(theme_bw())
#lreg.gglot <- abline(lm(raw.data$mix_1_Ct ~ raw.data$mix_2_Ct))

plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), level = 0.95, size = 0.5, color = "#4c4c4c") +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) +
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2a") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
#confidenceEllipse(lm(raw.data$mix_1_Ct ~ raw.data$mix_2_Ct), add = T, levels = 0.95)
#geom_line(lreg.ggplot)

ggsave(file="2ctB.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

#Same plot with differnt color on confidence ellipse
plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), level = 0.95, size = 1, color = "#b2b2b2") +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) +
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2b") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
#confidenceEllipse(lm(raw.data$mix_1_Ct ~ raw.data$mix_2_Ct), add = T, levels = 0.95)
#geom_line(lreg.ggplot)

ggsave(file="2ctB2.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

#Same plot with differnt color on confidence ellipse
plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), geom = "polygon", level = 0.95, size = 0.5, color = "#4c4c4c", alpha=0.2) +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) +
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2c") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))

ggsave(file="2ctB3.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

#Same plot with differnt color on confidence ellipse
plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.2) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), geom = "polygon", level = 0.95, size = 0.5, color = "#4c4c4c", alpha=0.2) +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) +
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2d") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))

ggsave(file="2ctB4.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

#Same plot with differnt color on confidence 
custom.grey <- c("#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2")

plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.2) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), geom = "polygon", level = 0.95, size = 0.2, color = "#4c4c4c", alpha=0.25) + scale_fill_manual(aesthetics = "fill", values = c("#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2")) +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) + 
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2e") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))

plot.qpcr$layers = rev(plot.qpcr$layers)

ggsave(file="2ctB5.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

Adjusting layering of the ellipses

#Same plot with differnt color on confidence 

df = raw.data
df = df[complete.cases(df),]
custom.grey <- c("#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2")

plot.qpcr <- ggplot(df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.2) + 
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), geom = "polygon", level = 0.95, size = 0, color = "#4c4c4c", alpha=0.5) + scale_fill_manual(aesthetics = "fill", values = c("#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2","#b1b2b2")) +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) + 
  labs(y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", title = "figure 2f") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))

plot.qpcr$layers = rev(plot.qpcr$layers)

ggsave(file="2ctB6.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

running the linear regression model to get equation

The results also contain the R-squared value.

model1 <- lm(df$mix_2_Ct ~ df$mix_1_Ct, na.action = "na.omit")
model1

## 
## Call:
## lm(formula = df$mix_2_Ct ~ df$mix_1_Ct, na.action = "na.omit")
## 
## Coefficients:
## (Intercept)  df$mix_1_Ct  
##     0.02848      0.99621

summary.lm(model1)

## 
## Call:
## lm(formula = df$mix_2_Ct ~ df$mix_1_Ct, na.action = "na.omit")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9266 -0.2666 -0.0021  0.2635  3.6047 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.028482   0.037779   0.754    0.451    
## df$mix_1_Ct 0.996209   0.001908 522.199   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5934 on 2212 degrees of freedom
## Multiple R-squared:  0.992,  Adjusted R-squared:  0.9919 
## F-statistic: 2.727e+05 on 1 and 2212 DF,  p-value: < 2.2e-16

paste("Ct-mix2 = 0.02848 + 0.99621*Ct-mix1")

## [1] "Ct-mix2 = 0.02848 + 0.99621*Ct-mix1"

Figure 3: plotting the prediction interval

temp_var <- predict(model1, interval="prediction")
new_df <- cbind(df, temp_var)

plot.pred <- ggplot(new_df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  geom_line(aes(y=lwr), color = "red", linetype = "dashed", fullrange = T) +
  geom_line(aes(y=upr), color = "red", linetype = "dashed") +
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct, fill = copies.33nl.well), level = 0.95, size = 0.5, color = "#4c4c4c") +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) + labs( 
       y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1") + 
  guides(color = guide_legend(override.aes = list(size=5), title="copies per well")) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
ggsave(file="linear.regression.confidence.intervals.pdf", device="pdf",plot=plot.pred, width=8, height=6)
plot(plot.pred)

temp_var <- predict(model1, interval="prediction")
new_df <- cbind(df, temp_var)

plot.pred <- ggplot(new_df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  geom_line(aes(y=lwr), color = "red", linetype = "dashed", fullrange = T) +
  geom_line(aes(y=upr), color = "red", linetype = "dashed") +
  stat_ellipse(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), level = 0.95, size = 0.5, color = "#4c4c4c") +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) + labs(subtitle="subtitle", 
       y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", 
       title="figure 3d") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
plot(plot.pred)

temp_var <- predict(model1, interval="prediction")
new_df <- cbind(df, temp_var)

plot.pred <- ggplot(new_df, aes(x=mix_1_Ct, y=mix_2_Ct, color = copies.33nl.well)) + xlim(0,40) + ylim(0,40) +
  geom_point(size = 0.5) + 
  geom_line(aes(y=lwr), color = "red", linetype = "dashed", fullrange = T) +
  geom_line(aes(y=upr), color = "red", linetype = "dashed") +
  geom_smooth(inherit.aes = F, aes(x=mix_1_Ct, y=mix_2_Ct), method = glm, se = TRUE, level = .99, size= 0.2, color = "red", fill = "blue", size = 1, fullrange = TRUE) + labs(subtitle="subtitle", 
       y="Ct RD-qChip Plate 2", 
       x="Ct RD-qChip Plate 1", 
       title="figure 3e") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"))
plot(plot.pred)

Box Plots

library(ggplot2)
library(scales)

df = raw.data
df = df[complete.cases(df),]

theme_set(theme_bw())
plot.qpcr <- ggplot(df, aes(x=copies.33nl.well, y=mix_1_Ct)) +
  geom_jitter(size = 0.2, aes(color = copies.33nl.well), position = position_jitter(width=c(.5,0.2))) +
  stat_boxplot(geom ="errorbar", width = 0.25) + 
  geom_boxplot(outlier.alpha = 0, alpha = .1, width = 0.5) +
  stat_summary(fun.y=mean, colour="black", geom="point", 
               shape=21, fill = "white", size=2,show_guide = FALSE) +
  labs(y="Ct plate 1", 
       x="copies per well") + guides(color = guide_legend(override.aes = list(size=5))) + 
  labs(color='copies per well') +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"),axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file="4ctA.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr)

library(ggplot2)
library(scales)

df = raw.data
df = df[complete.cases(df),]

#with means
theme_set(theme_bw())
plot.qpcr1 <- ggplot(df, aes(x=copies.33nl.well, y=mix_1_Ct)) +
  geom_jitter(size = 0.1, aes(color = copies.33nl.well)) +
  stat_boxplot(geom ="errorbar", width = 0.25) + 
  geom_boxplot(outlier.alpha = 0, alpha = 0) +
  ylim(0,35) +
    stat_summary(fun.y=mean, colour="black", geom="point", 
               shape=21, fill = "white", size=2,show_guide = FALSE) +
  labs(y="Ct plate 1", 
       x="copies/33nl hole", title = "figure 4f") + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"),axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file="4ctF.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr1)

theme_set(theme_bw())
plot.qpcr2 <- ggplot(df, aes(x=copies.33nl.well, y=mix_2_Ct)) +
  geom_jitter(size = 0.1, aes(color = copies.33nl.well)) +
  stat_boxplot(geom ="errorbar", width = 0.25) + 
  geom_boxplot(outlier.alpha = 0, alpha = 0) +
  ylim(0,35) +
    stat_summary(fun.y=mean, colour="black", geom="point", 
               shape=21, fill = "white", size=2,show_guide = FALSE)+
  labs(y="Ct plate 2", 
       x="copies per 33nl well", title = "figure 4g") + guides(color = guide_legend(override.aes = list(size=5))) + 
  labs(color='copies/33nl hole') +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"),axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file="4ctG.pdf", device="pdf",plot=plot.qpcr, width=8, height=6)
plot(plot.qpcr2)

#slightly different color scheme
plot.qpcr1 <- ggplot(df, aes(x=copies.33nl.well, y=mix_1_Ct)) +
  geom_jitter(size = 0.2, aes(color = copies.33nl.well)) +
  stat_boxplot(geom ="errorbar", width = 0.25) + 
  geom_boxplot(outlier.alpha = 0, alpha = 0) +
  ylim(0,35) +
    stat_summary(fun.y=mean, colour="black", geom="point", 
               shape=21, fill = "white", size=2,show_guide = FALSE) +
  labs(y="Ct plate 1", 
       x="copies per well", title = NULL) + guides(color = guide_legend(override.aes = list(size=5))) +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"),axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file="4ctH.pdf", device="pdf",plot=plot.qpcr1, width=8, height=6)
plot(plot.qpcr1)

theme_set(theme_bw())
plot.qpcr2 <- ggplot(df, aes(x=copies.33nl.well, y=mix_2_Ct)) +
  geom_jitter(size = 0.2, aes(color = copies.33nl.well)) +
  stat_boxplot(geom ="errorbar", width = 0.25) + 
  geom_boxplot(outlier.alpha = 0, alpha = 0) +
  ylim(0,35) +
    stat_summary(fun.y=mean, colour="black", geom="point", 
               shape=21, fill = "white", size=2,show_guide = FALSE)+
  labs(y="Ct plate 2", 
       x="copies per 33nl well", title = NULL) + guides(color = guide_legend(override.aes = list(size=5))) + 
  labs(color='copies per well') +
  theme(axis.text=element_text(size=12,face="bold"),
        axis.title=element_text(size=14,face="bold"),axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file="4ctI.pdf", device="pdf",plot=plot.qpcr2, width=8, height=6)
plot(plot.qpcr2)