Examples/PlotDistribution.R

In matsengrp/sumrep: Summary statistics for B cell receptor (BCR) repertoires

library(viridis)

plotDistribution <- function(function_string,
                             datasets,
                             column_name,
                             dat_names,
                             group_names,
                             filename=paste0(function_string, ".pdf")) {

    summary_function <- get(function_string)
    d_full <- {}
    num_datasets <- length(datasets)
    for(i in 1:num_datasets) {
        density_i <- datasets[[i]][[column_name]] %>%
            summary_function
        d_full <- rbind(d_full,
                        data.table(Value=density_i,
                                   Dataset=dat_names[[i]],
                                   Group=group_names[[i]])
                        )
    }

    ggplot(d_full, aes(x=Value, colour=Dataset, lty=Group)) + 
        geom_density() +
        scale_color_viridis(discrete=TRUE) +
        ggtitle(paste("Densities of", function_string))

    ggsave(filename, width=10, height=6)
}

dats <- list(p_f1$annotations,
             p_f2$annotations,
             p_g1$annotations,
             p_f1_sim$annotations,
             p_f2_sim$annotations,
             p_g1_sim$annotations)
groups <- c(rep("Observed", 3), rep("Simulated", 3))

plotDistribution(function_string="getGCContentDistribution",
                 datasets=dats,
                 column_name="sequence",
                 dat_names=c("Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1",
                             "Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1"),
                 group_names=groups
                )

                 
plotDistribution(function_string="getHydrophobicityDistribution",
                 datasets=dats,
                 column_name="junction",
                 dat_names=c("Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1",
                             "Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1"),
                 group_names=groups
                )

plotDistribution(function_string="getGRAVYDistribution",
                 datasets=dats,
                 column_name="junction",
                 dat_names=c("Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1",
                             "Individual 1, time 1",
                             "Individual 1, time 2",
                             "Individual 2, time 1"),
                 group_names=groups
                )