Analysis

Read dataset and preprocess.

Remove mAb 135 because it is missing Neut_micro.

df <- read.csv("./data/master_log.csv", na.strings=c("?", "nd"), strip.white=TRUE, row.names=c("Ab"))
df <- df[!rownames(df) %in% c("135"),]

Convert Epitope Class to a factor.

df$Epitope_Class <- factor(df$Epitope_Class)

Drop columns with missing values.

df <- subset(df, select=-c(Protect_binary, Epitope_Class_ELISA, Endotoxin, Endotoxin.1, mW_Loss, aTTD, Epitope.Class..assigning.method.)) #Remove columns with empty values.
df[,"Escape..code"] <- as.factor(df[,"Escape..code"])
df[,"Makona.binding"] <- as.factor(df[,"Makona.binding"])

Make unNeutFrac, Neut_dVP30 and Neut_VSV correlate positively with Protection.

df[,"unNeutFrac"] <- 1 - df[,"unNeutFrac"]
df[,"Neut_dVP30"] <- 1- df[,"Neut_dVP30"]
df[,"Neut_VSV"] <- (100 - df[,"Neut_VSV"])/100

Separate Human IgG1 and Mouse IgG1

levels(df[,"Isotype"]) <- c(levels(df[,"Isotype"]), "HumanIgG1")
df[df[,"Isotype"]=="IgG1" & df[,"Species"]=="human", "Isotype"] <- "HumanIgG1"
levels(df[,"Isotype"]) <- c(levels(df[,"Isotype"]), "MouseIgG1")
df[df[,"Isotype"]=="IgG1" & df[,"Species"]=="mouse", "Isotype"] <- "MouseIgG1"
df[,"Isotype"] <- droplevels(df[,"Isotype"])

Convert polyfunctionality to numeric type.

df[,"Polyfunctionality"] <- as.numeric(df[,"Polyfunctionality"])

Add Epitope Tier columns

df[df$Epitope_Class %in% c("Cap", "GP1/Head", "Mucin"), "Epitope_Tier"] <- "Tier1"
df[df$Epitope_Class %in% c("Base", "GP1/Core", "Fusion"), "Epitope_Tier"] <- "Tier2"
df[df$Epitope_Class %in% c("GP1/2", "HR2"), "Epitope_Tier"] <- "Tier3"
df[df$Epitope_Class %in% c("Unknown"), "Epitope_Tier"] <- "TierUnknown"

Create mAB id columns

df[,"id"] <- rownames(df)
head(df)
##   Protection Epitope_Class  Neut_VSV unNeutFrac Neut_micro Neut_dVP30
## 1        0.5           Cap 0.0000000        0.0       0.00          0
## 2        0.0         Mucin 0.0000000        0.0       0.00          0
## 3        0.1           Cap 0.0000000        0.0       0.20          0
## 4        0.2         Mucin 0.0000000        0.0       0.00          0
## 5        0.3           Cap 0.9875636        0.7       0.62          0
## 6        0.4         Mucin 0.0000000        0.0       0.00          0
##   GP_EC50 sGP_EC50 Species   Isotype Round Cross.reactivity
## 1    1.38     0.61   human HumanIgG1    R4             None
## 2    0.91  1000.00   human HumanIgG1    R6             None
## 3  175.80    22.19   mouse     IgG2a    R2             None
## 4    0.68  1000.00   human HumanIgG1    R2             None
## 5   44.24     7.45   mouse     IgG2a    R1             None
## 6    4.31  1000.00   mouse     IgG2a    R2             None
##   Polyfunctionality huADCP mADCP CD107a  IFNg MIP.1b huADNP mADNP
## 1                 6  43.32 73.93  33.39 32.10  56.48 154.45 40.68
## 2                 6  25.58 66.50  40.54 39.10  74.28 126.64 35.44
## 3                 2  38.01 45.89   8.98  5.95  15.10  16.45 39.01
## 4                 6  39.95 62.92  45.94 41.15  84.93 101.97 31.96
## 5                 2  51.91 44.98   9.53  6.37  16.39  19.28 32.63
## 6                 3  31.50 58.63  11.15  6.04  18.59  42.10 41.57
##   sGP_huADCP sGP_huADNP Total_G0 Total_G1 Total_G2 Total_Fucose
## 1     274.70      57.57    85.06     4.34    10.60        12.86
## 2     162.02      46.01    72.10    17.31    10.59        56.49
## 3     269.58      27.54    55.52    31.83    12.66        98.14
## 4     164.10      31.82    88.83     2.71     8.47        12.47
## 5     214.98      31.26    66.79    18.57    14.64        97.03
## 6     107.00      25.73    61.53    25.17    13.30        98.24
##   Total_Bisecting Total_mono_SA Total_SA    G0   G0F G0FB   G1  G1.  G1F
## 1            2.51          8.44     8.44 83.41  1.62 0.03 0.12 1.24 0.02
## 2            1.31         11.21    11.21 40.39 31.35 0.36 0.43 0.61 2.53
## 3            0.81         13.02    13.02  0.67 54.38 0.47 0.33 0.80 2.46
## 4            1.85          6.78     6.78 83.24  5.18 0.41 0.17 1.29 0.15
## 5            2.05         17.02    17.02  0.21 66.12 0.47 0.45 1.05 1.43
## 6            0.77         13.53    13.53  0.26 60.83 0.44 0.37 0.45 2.06
##    G1F. G1S1F G1FB  G2F G2FB G2S1 G2S1F G2S1B G2S1FB      Escape..code
## 1  0.01  2.73 0.23 4.80 0.08 1.10  2.44  1.27   0.90 No Escape Mutants
## 2 10.23  3.26 0.26 2.58 0.06 1.64  5.69  0.44   0.18 No Escape Mutants
## 3 26.51  1.62 0.11 1.20 0.06 0.05 11.17  0.01   0.16 No Escape Mutants
## 4  0.21  0.84 0.06 2.35 0.18 2.03  2.69  0.81   0.41 No Escape Mutants
## 5 11.76  3.75 0.13 1.00 0.37 0.23 11.96  1.03   0.05 No Escape Mutants
## 6 20.82  1.45 0.03 1.19 0.03 0.42 11.39  0.27   0.00 No Escape Mutants
##   Makona.binding sGP.binder Epitope_Tier id
## 1       Not done       TRUE        Tier1  1
## 2       Not done      FALSE        Tier1  2
## 3       Not done       TRUE        Tier1  3
## 4       Not done      FALSE        Tier1  4
## 5       Not done       TRUE        Tier1  5
## 6       Not done      FALSE        Tier1  6

Univariate Analysis

Calculate Spearman’s rank correlation coefficient between “Protection” and other columns. Plot only statistically significant correlations(significance value = 0.05).

library(ggplot2)

## Correlation Strength
co.pvalue <- sapply(colnames(df), function(x){
    if(is.numeric(df[,x])){
        t <- cor.test(df[,"Protection"], df[,x], method="spearman", exact=FALSE)
        t$p.value
    }
})
co.estimate <- sapply(colnames(df), function(x){
    if(is.numeric(df[,x])){
        t <- cor.test(df[,"Protection"], df[,x], method="spearman", exact=FALSE)
        as.numeric(t$estimate)
    }
})

co <- data.frame(unlist(co.pvalue), unlist(co.estimate))
colnames(co) <- c("pval", "corr")
co$pval <- p.adjust(co$pval, method="BH")
co <- co[co$pval <= 0.05,]

rownames(co) <- sapply(rownames(co), function(x){
    gsub(".rho", "", x)
})
co[,"X1"] <- "Protection"
co[,"X2"] <- rownames(co)
co <- co[with(co, order(-abs(corr))),]
co <- co[co[,"X2"]!="Protection",]

ggplot(co, aes(reorder(X2, abs(corr)), abs(corr))) + geom_bar(stat="identity", aes(fill=corr > 0)) + scale_fill_manual(name="Spearman's Rank Correlation", labels = c("FALSE" = "Negative", "TRUE" = "Positive"), values = c('#ff4c4c', '#4c4cff')) + xlab("") + ylab("Correlation Strength") + theme_bw()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ coord_flip()

Identify differences in features grouped by Epitope_Class

Kruskal-Wallis test for numeric features. Chi squared test for categorical features. Significance Level = 0.05

t <- df

pval <- c()
col <- c()
fpval <- c()
fcol <- c()
for(i in colnames(t)){
    if(class(t[,i])=="numeric" & i!="Epitope_Class"){ #For numeric features perform Kruskal–Wallis test.
        tdf <- data.frame(feature=t[,i], group= t[,"Epitope_Class"])
        f <- kruskal.test(feature~group, data=tdf)
        pval <- c(pval, f$p.value)
        col <- c(col, i)
    } else if(class(t[,i])!="numeric" & i!="Epitope_Class") { #For categorical features perform Chi Squared test.
        t[,i] <- as.factor(t[,i])
        print(i)
        set.seed(1)
        p <- chisq.test(as.factor(t[,i]), t[,"Epitope_Class"], simulate.p.value=TRUE)
        fpval <- c(fpval, p$p.value)
        fcol <- c(fcol, i)
    }
}
## [1] "Species"
## [1] "Isotype"
## [1] "Round"
## [1] "Cross.reactivity"
## [1] "Escape..code"
## [1] "Makona.binding"
## [1] "sGP.binder"
## [1] "Epitope_Tier"
## [1] "id"

Multiple Hypothesis Correction by false discovery rate(Benjamani-Hochberg Method). Significance Level = 0.05

corpval <- p.adjust(unlist(pval), method="fdr")
fcorpval <- p.adjust(unlist(fpval), method="fdr")

epdiff <- data.frame(col, pval, corpval)
fepdiff <- data.frame(col=fcol, pval=fpval, corpval=fcorpval)

epdiff[epdiff[,"corpval"]<0.05,]
##                  col         pval      corpval
## 1         Protection 5.923109e-08 3.300018e-07
## 2           Neut_VSV 7.478246e-15 7.291289e-14
## 3         unNeutFrac 5.360700e-15 6.968910e-14
## 4         Neut_micro 4.077704e-12 2.650508e-11
## 5         Neut_dVP30 1.901841e-12 1.483436e-11
## 6            GP_EC50 2.068507e-03 5.762271e-03
## 7           sGP_EC50 6.359297e-26 2.480126e-24
## 8  Polyfunctionality 4.230926e-03 1.100041e-02
## 9             huADCP 4.434503e-06 1.921618e-05
## 10             mADCP 6.597975e-05 2.573210e-04
## 14            huADNP 1.132544e-03 3.397632e-03
## 15             mADNP 2.804426e-04 9.942967e-04
## 16        sGP_huADCP 9.127136e-17 1.779791e-15
## 17        sGP_huADNP 3.970938e-07 1.935833e-06
## 19          Total_G1 1.906499e-02 4.130749e-02
## 26               G0F 1.745935e-02 4.005380e-02
## 30               G1F 2.372277e-02 4.869411e-02
## 31              G1F. 5.059088e-03 1.233153e-02
## 33              G1FB 6.877308e-04 2.235125e-03
fepdiff[fepdiff[,"corpval"]<0.05,]
##                col         pval     corpval
## 1          Species 0.0014992504 0.003373313
## 2          Isotype 0.0039980010 0.005140287
## 3            Round 0.0004997501 0.001499250
## 4 Cross.reactivity 0.0039980010 0.005140287
## 5     Escape..code 0.0069965017 0.007871064
## 6   Makona.binding 0.0024987506 0.004497751
## 7       sGP.binder 0.0004997501 0.001499250
## 8     Epitope_Tier 0.0004997501 0.001499250

Use Dunn’s test for pairwise comparisons between groups with statistically significant differences Create box plots of groups annotated with statistically significant differences, labelled outliers and mean(red dot). Pretty plot groups without annotations and using color_scheme without annotations for figures.

library(dunn.test)
library(ggsignif)
library(dplyr)
## devtools::install_github("thomasp85/patchwork")
library(patchwork)

color_scheme_epitope_class <- c('#e03a3e','#0B9CD8','#963e67','#F58220','#61bb46','#6D6E71','#fdb724','#2e3192','#000000')
names(color_scheme_epitope_class) <- c('Base','Cap','Fusion','GP1/Core','GP1/Head','GP1/2','HR2','Mucin','Unknown')

generate_annotations <- function(t, ymax, diff){
    start = c()
    end = c()
    y = c()
    label=c()
    for(i in rownames(t)){
        if(t[i ,"P.adjusted"] <= 0.05){
            ymax = ymax + (diff/length(names(t)))
            v <- unlist(strsplit(as.character(t[i, "comparisons"]), " - "))
            start <- c(start, v[1])
            end <- c(end, v[2])
            y <- c(y, ymax)
            label <- c(label, signif(as.numeric(t[i, "P.adjusted"]), 3))
        }
    }
    return(data.frame(start, end, y, label))
}

is_outlier <- function(x) {
    return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

fun_mean <- function(x){
    return(data.frame(y=mean(x),label=signif(mean(x,na.rm=T), 2)))
}

t.f <- df

cols <- epdiff[epdiff[,"corpval"]<0.05,"col"]
plts <- list()

orderEpitopeClass <- c("Unknown", "Mucin", "GP1/2", "GP1/Core", "Cap", "HR2","Fusion","Base","GP1/Head")

## Dunn's test
for(i in cols){
    t.f[,"Epitope_Class"] <- factor(t.f$Epitope_Class, levels=orderEpitopeClass)
    temp <- dunn.test(t.f[,i], t.f[,"Epitope_Class"], method="bh")
    temp <- as.data.frame(temp)
    temp <- temp[with(temp, order(comparisons)),]
    anndf <- generate_annotations(temp, max(t.f[,i]), (range(t.f[,i])[2] - range(t.f[,i])[1]))
    outliers <- t.f %>% group_by(Epitope_Class) %>% mutate(outlier = ifelse(is_outlier(eval(as.name(paste(i)))), id, ""))
    t.f[,"outlier_label"] <- outliers[["outlier"]]
    t.f[,paste("outlier_label", i, sep="_")] <- outliers[["outlier"]]
    p <- ggplot(t.f, aes_string(x="Epitope_Class", y=`i`)) + geom_boxplot() + stat_summary(fun.y = mean, geom = "point", size = 1, color="red", fill="red", shape=23) + stat_summary(fun.data = fun_mean, geom="text", vjust=1.5, color="red") + geom_text(aes(label=outlier_label), hjust=-0.3)+theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_signif(data=anndf, aes(xmin=start, xmax=end, annotations=label, y_position=y), textsize = 3, vjust = -0.2, manual=TRUE)
    plts[[length(plts)+1]] <- p
    medians <- aggregate(t.f[,i], by=list(t.f[,"Epitope_Class"]), FUN=median)
    colnames(medians) <- c("Group", "x")
    medians[,"Group"] <- factor(medians$Group,levels=orderEpitopeClass)
    p <- ggplot(t.f, aes_string(x="Epitope_Class", y=`i`))+geom_jitter(aes(colour=Epitope_Class), stroke=0, position=position_jitter(width = 0.1, height=0), size=4, alpha=0.7)+scale_colour_manual(guide=FALSE,name = "Epitope Class",values = color_scheme_epitope_class) + stat_summary(fun.data="median_hilow", geom="linerange", size=0.2, colour="black")  + theme_bw()+ theme(panel.border = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line( size=.3, color="#f5f5f5"), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.text.x = element_text(angle = 45, hjust = 1)) +xlab("Epitope Class")+ geom_segment(aes(x=as.numeric(Group)-0.3, xend=as.numeric(Group)+0.3, y=x, yend=x), data=medians)+xlab("")+ylab(i)+ scale_y_continuous(labels = scales::percent)
    plts[[length(plts)+1]] <- p
}

p <- plts[[1]]
for(x in plts[-1]){
    p <- p+ x
}
p + plot_layout(ncol = 2)

Identify strong, partial and weak mAbs based on neutralization readouts

We use k-means clustering of mABs using neutralization readouts - Neut_micro, Neut_VSV, Neut_dVP30 and unNeutFrac. To specify the number of clusters, ‘k’, we use the elbow method.

df.neut <- df[,c("Neut_micro", "Neut_VSV", "Neut_dVP30", "unNeutFrac")]
df.kmeans <- data.frame(clusters=c(1:10), wss=rep(c(0),times=10))
for (i in 1:10){
    set.seed(112358)
    df.kmeans[i, "wss"] <- sum(kmeans(df.neut, centers=i)$withinss)
}

ggplot(data=df.kmeans, aes(x=clusters, y=wss, group=1)) +
  geom_line()+
  geom_point()

Let’s select k = 5 based and use t-SNE to visualize the 5 clusters.

set.seed(112358)
df.kmeans <-  kmeans(df.neut, centers=5)
df.nrneut <- df.neut

library(Rtsne)
library(RColorBrewer)
set.seed(112358)
df.tsne <- Rtsne(as.matrix(df.nrneut), check_duplicates = FALSE, max_iter = 5000, perplexity = 20)
df.nrneut[,"x"] <- df.tsne$Y[,1]
df.nrneut[,"y"] <- df.tsne$Y[,2]
df.nrneut[,"cluster"] <- df.kmeans$cluster
df.nrneut$cluster <- as.factor(df.nrneut[,"cluster"])
ggplot(data=df.nrneut, aes(x= x, y = y, color=cluster)) + geom_point()+ scale_color_brewer("Cluster", palette="Dark2")

Let’s now plot the four neutralization readouts grouped by clusters.

library(reshape2)
df.nrneut$id <- factor(df$id)
df.nrneut$epclass <- df[,"Epitope_Class"]
df.nrneut[,"Protection"] <- df[,"Protection"]
df.nrneut$cluster <- factor(df.nrneut$cluster, levels=c(3,4,5,1,2))
## df.nrneut$id <- as.numeric(df.nrneut$id)
df.nrneut <- df.nrneut[with(df.nrneut, order(cluster, -id)),]
## Warning in Ops.factor(id): '-' not meaningful for factors
df.nrneut$id <- factor(df.nrneut$id, levels=df.nrneut$id)
df.neut <- melt(subset(df.nrneut, select=-c(x,y)))
## Using cluster, id, epclass as id variables
df.neut$id <- factor(df.neut$id, levels=df.neut$id)
## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated
df.neut$xlabel <- paste(df.neut, df.neut$id)
p_neut <- ggplot() +geom_point(data=df.neut, size=3, aes(x=id, y=variable)) + geom_point(data=df.neut, size=2.5, aes(colour=value, x=id, y=variable)) +scale_colour_gradientn(colors=brewer.pal(9,"YlOrRd")) + theme(panel.background = element_rect(fill = '#FFFFFF', colour = '#FFFFFF'), axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5, size=8), axis.text.y = element_text(size=14), text=element_text(size=12))+ xlab("Antibody") + ylab("")
y1 <- c()
y2 <- c()
x1 <- c()
x2 <- c()
ep <- c()
for(i in levels(df.neut[,"cluster"])){
    u <- unique(df.neut[,"id"][df.neut["cluster"]==i])
    x1 <- c(x1, as.numeric(u[1])-0.45)
    x2 <- c(x2, as.numeric(tail(u,n = 1))+0.45)
    y1 <- c(y1, min(as.numeric(df.neut$variable))-0.5)
    y2 <- c(y2, max(as.numeric(df.neut$variable))+0.5)
    ep <- c(ep, i)
}
## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated
## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated
shade <- data.frame(x1, x2, y1, y2, ep)
p_neut <- p_neut+geom_rect(data=shade,mapping = aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2), color=brewer.pal(length(shade[,"ep"]), "Set1"), alpha=0, fill="transparent", size=0.5)

## Add Annotations
for(i in rownames(shade)){
    p_neut <- p_neut + annotate("text",x=shade[i, "x1"],y=shade[i, "y2"], angle=0,hjust=-0.5,vjust=1.2,label=shade[i, "ep"],size=4)
}
p_neut + coord_flip()
## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

## Warning in `levels<-`(`*tmp*`, value = if (nl == nL) as.character(labels)
## else paste0(labels, : duplicated levels in factors are deprecated

Based on the plot above, let’s classify the clusters as following,

  • 2 -> “Weak”
  • 1,4,5 -> “Partial”
  • 3 -> “Strong”
df[,"cluster"] <- df.nrneut$cluster
df[df[,"cluster"]=="2","groups"] <- "weak"
df[df[,"cluster"]=="1","groups"] <- "partial"
df[df[,"cluster"]=="4","groups"] <- "partial"
df[df[,"cluster"]=="5","groups"] <- "partial"
df[df[,"cluster"]=="3","groups"] <- "strong"

Let’s now use Dunn’s test to check if the difference in means of the neutralization readouts between all three groups is statistically significant, thus validating our previous step.

t.f <- df
plts <- list()

groups_color_scheme <- c("#D75466", "#ffc107", "#17a2b8")
cols <- c("Neut_micro", "unNeutFrac", "Neut_dVP30", "Neut_VSV", "Protection")
orderGroups <- c("weak", "partial", "strong")
for(i in cols){
    t.f[,"groups"] <- factor(t.f[,"groups"], levels = orderGroups)
    temp <- dunn.test(t.f[,i], t.f[,"groups"], method="bh")
    temp <- as.data.frame(temp)
    temp <- temp[with(temp, order(comparisons)),]
    anndf <- generate_annotations(temp, max(t.f[,i]), (range(t.f[,i])[2] - range(t.f[,i])[1]))
    ## Boxplot
    outliers <- t.f %>% group_by(groups) %>% mutate(outlier = ifelse(is_outlier(eval(as.name(paste(i)))), id, ""))
    t.f[,"outlier_label"] <- outliers[["outlier"]]
    t.f[,paste("outlier_label", i, sep="_")] <- outliers[["outlier"]]
    p <- ggplot(t.f, aes_string(x="groups", y=`i`)) + geom_boxplot() + stat_summary(fun.y = mean, geom = "point", size = 1, color="red", fill="red", shape=23) + stat_summary(fun.data = fun_mean, geom="text", vjust=1.5, color="red") + geom_text(aes(label=outlier_label), hjust=-0.3)+theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_signif(data=anndf, aes(xmin=start, xmax=end, annotations=label, y_position=y), textsize = 3, vjust = -0.2, manual=TRUE)
    plts[[length(plts)+1]] <- p
    medians <- aggregate(t.f[,i], by=list(t.f[,"groups"]), FUN=median)
    colnames(medians) <- c("Group", "x")
    medians[,"Group"] <- factor(medians$Group,levels=orderGroups)
    p <- ggplot(t.f, aes_string(x="groups", y=`i`))+geom_jitter(aes(colour=groups), stroke=0, position=position_jitter(width = 0.1, height=0), size=4, alpha=0.7)+scale_colour_manual(guide=FALSE,name = "Groups",values = groups_color_scheme) + stat_summary(fun.data="median_hilow", geom="linerange", size=0.2, colour="black")  + theme_bw()+ theme(text = element_text(size=14), axis.text = element_text(size=14), panel.border = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line( size=.3, color="#f5f5f5"), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.text.x = element_text(angle = 45, hjust = 1)) +xlab("Groups")+ geom_segment(aes(x=as.numeric(Group)-0.3, xend=as.numeric(Group)+0.3, y=x, yend=x), data=medians)+xlab("Groups")+ylab(i)
    plts[[length(plts)+1]] <- p
}
p <- plts[[1]]
for(x in plts[-1]){
    p <- p+ x
}
p + plot_layout(ncol = 2)

We see that at a significance level of 0.05, the difference in the means of four neutralization readouts across the three groups is statistically significant thus validating out previous step of classifying the clusters int ostrong, weak and partial neutralizers.