Diversity report for diversity outbred mice

After finishing 06_final_pr_apr_69K.R, 07_do_diversity_report.R, all the output will be used plot DO Diversity Report for 12 batches of DO mice

library

# Load packages
library(qtl2)
library(table1)
library(tidyr)
library(dplyr)
library(data.table)
library(foreach)
library(doParallel)
library(parallel)
library(abind)
library(gap)
library(regress)
library(lme4)
library(abind)
library(ggplot2)
library(vcd)
library(MASS)
library(plotly)
library(colorspace)
options(stringsAsFactors = FALSE)
source("code/reconst_utils.R")

#Summary

load("data/Jackson_Lab_12_batches/gm_DO3173_qc.RData")#gm_after_qc

# make dataset with a few variables to summarize
table1 <- gm_after_qc$covar %>% 
  select(Name = name, 
         Sex  = sex, 
         Generation = ngen) %>%
  mutate(Sex = case_when(
    Sex == "F" ~ "Female",
    Sex == "M" ~ "Male"
  ))

# summarize the data 
table1(~ Generation | Sex, data=table1)

#Founder contributions

load("data/Jackson_Lab_12_batches/fp_DO3173.RData") #fp and fp_summary object
#change order of level in gen
fp$gen <- factor(fp$gen,levels = c(21,22,23,25,29,30,31,32,33,34,35,36))

#summarize per generation per chromosome
fp_summary = fp %>% group_by(chr, founder, gen) %>%
  summarize(mean = round(100*mean(prop), 2),
            sd   = round(100*sd(prop), 2))

`summarise()` regrouping output by 'chr', 'founder' (override with `.groups` argument)

#Stackbar plot
#summarize per chromosome across generation
pdf(file = "data/Jackson_Lab_12_batches/stackbar_mean_prop_across_all_gen.pdf",width = 16)
p01 <- fp %>% group_by(chr, founder) %>%
  summarise(grand_mean = round(100*mean(prop), 2)) %>%
  ggplot(aes(x = chr, y = grand_mean, fill = founder)) +
    geom_bar(stat="identity",
             width=1) +
    geom_text(aes(label = paste0(grand_mean)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    ylab("Mean percentage across generations") +
    xlab("Chromosome") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())

`summarise()` regrouping output by 'chr' (override with `.groups` argument)

p01
dev.off()

png 
  2 

p01

Past versions of “Founder contributions-1.png”

#Stackbar plot
#summarize per chromosome across generation
pdf(file = "data/Jackson_Lab_12_batches/stackbar_mean_prop_across_all_chr.pdf",width = 16)
p02 <- fp %>% group_by(gen, founder) %>%
  summarise(grand_mean = round(100*mean(prop), 2), 
            grand_sd   = round(100*sd(prop), 2)) %>%
  ggplot(aes(x = gen, y = grand_mean, fill = founder)) +
    geom_bar(stat="identity",
             width=0.99) +
    geom_text(aes(label = paste0(grand_mean, " ± ", grand_sd)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    ylab("Mean percentage across all chromosomes") +
    xlab("Generation") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())

`summarise()` regrouping output by 'gen' (override with `.groups` argument)

p02
dev.off()

png 
  2 

p02

Past versions of “Founder contributions-2.png”

#stackbar_prop_across_gen
for(c in c(1:19, "X")){
  #print(c)
  p <- ggplot(data = fp_summary[fp_summary$chr == c,], aes(x = gen, y = mean, fill = founder)) +
    geom_bar(stat="identity",
             width=1) +
    geom_text(aes(label = paste0(mean," ± ", sd)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    labs(title = paste0("Chr ", c)) +
    ylab("Percentage") +
    xlab("Generation") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())
  print(p)
}

Past versions of “Founder contributions-3.png”

Past versions of “Founder contributions-4.png”

Past versions of “Founder contributions-5.png”

Past versions of “Founder contributions-6.png”

Past versions of “Founder contributions-7.png”

Past versions of “Founder contributions-8.png”

Past versions of “Founder contributions-9.png”

Past versions of “Founder contributions-10.png”

Past versions of “Founder contributions-11.png”

Past versions of “Founder contributions-12.png”

Past versions of “Founder contributions-13.png”

Past versions of “Founder contributions-14.png”

Past versions of “Founder contributions-15.png”

Past versions of “Founder contributions-16.png”

Past versions of “Founder contributions-17.png”

Past versions of “Founder contributions-18.png”

Past versions of “Founder contributions-19.png”

Past versions of “Founder contributions-20.png”

Past versions of “Founder contributions-21.png”

Past versions of “Founder contributions-22.png”

pdf(file = "data/Jackson_Lab_12_batches/stackbar_prop_across_gen.pdf",width = 16)
for(c in c(1:19, "X")){
  #print(c)
  p <- ggplot(data = fp_summary[fp_summary$chr == c,], aes(x = gen, y = mean, fill = founder)) +
    geom_bar(stat="identity",
             width=1) +
    geom_text(aes(label = paste0(mean," ± ", sd)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    labs(title = paste0("Chr ", c)) +
    ylab("Percentage") +
    xlab("Generation") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())
  print(p)
}
dev.off()

png 
  2 

#stackbar_prop_across_chr
for(g in levels(fp_summary$gen)){
  #print(g)
  p <- ggplot(data = fp_summary[fp_summary$gen == g,], aes(x = chr, y = mean, fill = founder)) +
    geom_bar(stat="identity",
             width=1) +
    geom_text(aes(label = paste0(mean)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    labs(title = paste0("Generation ", g)) +
    ylab("Percentage") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.title.x = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())
  #print(p)
}

pdf(file = "data/Jackson_Lab_12_batches/stackbar_prop_across_chr.pdf", width = 12)
for(g in levels(fp_summary$gen)){
  #print(g)
  p <- ggplot(data = fp_summary[fp_summary$gen == g,], aes(x = chr, y = mean, fill = founder)) +
    geom_bar(stat="identity",
             width=1) +
    geom_text(aes(label = paste0(mean)), position = position_stack(vjust = 0.5)) +
    scale_fill_manual(values = CCcolors) +
    labs(title = paste0("Generation ", g)) +
    ylab("Percentage") +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.title.x = element_blank(),
          axis.line = element_line(colour = "black"),
          plot.title = element_text(hjust = 0.5),
          text = element_text(size=16),
          axis.title=element_text(size=16),
          legend.title=element_blank())
  print(p)
}
dev.off()

png 
  2 

#line plot
#plt <- htmltools::tagList()
for(c in unique(names(gm_after_qc$geno))){
  print(c)
  fp_subdata <- fp[fp$chr == c,]
  pp <- ggplot(data = fp_subdata,aes(pos, prop, group = gen, color = founder)) +
    geom_line(aes(linetype=gen)) +
    scale_linetype_manual(values=rep("solid",12)) +
    geom_hline(yintercept=0.125, linetype="dashed", color = "black", size = 0.25) +
    scale_color_manual(values = CCcolors) +
    facet_grid(founder~.) + 
    labs(title = paste0("Chr ", c)) + 
    theme(legend.position='none')
  print(pp)
  # Print an interactive plot
  # Add to list
  #plt[[c]] <- as_widget(ggplotly(pp, width = 1000, height = 1000))
}

[1] "1"

Past versions of “Founder contributions-23.png”

[1] "2"

Past versions of “Founder contributions-24.png”

[1] "3"

Past versions of “Founder contributions-25.png”

[1] "4"

Past versions of “Founder contributions-26.png”

[1] "5"

Past versions of “Founder contributions-27.png”

[1] "6"

Past versions of “Founder contributions-28.png”

[1] "7"

Past versions of “Founder contributions-29.png”

[1] "8"

Past versions of “Founder contributions-30.png”

[1] "9"

Past versions of “Founder contributions-31.png”

[1] "10"

Past versions of “Founder contributions-32.png”

[1] "11"

Past versions of “Founder contributions-33.png”

[1] "12"

Past versions of “Founder contributions-34.png”

[1] "13"

Past versions of “Founder contributions-35.png”

[1] "14"

Past versions of “Founder contributions-36.png”

[1] "15"

Past versions of “Founder contributions-37.png”

[1] "16"

Past versions of “Founder contributions-38.png”

[1] "17"

Past versions of “Founder contributions-39.png”

[1] "18"

Past versions of “Founder contributions-40.png”

[1] "19"

Past versions of “Founder contributions-41.png”

[1] "X"

Past versions of “Founder contributions-42.png”

#plt

#Average haplotype block size

load("data/Jackson_Lab_12_batches/recom_block_size.RData")
#Create an appropriately sized vector of names
nameVector <- unlist(mapply(function(x,y){ rep(y, length(x)) }, pos_ind_gen, names(pos_ind_gen)))
#Create the result
recom_block <- cbind.data.frame(unlist(pos_ind_gen), nameVector)
colnames(recom_block) <- c("sizeblock",
                           "ngen")
#remove 0
recom_block <- recom_block[recom_block$sizeblock != 0,]
recom_block$ngen <- factor(recom_block$ngen, levels = as.character(c(21:36)))
#mean
means <- aggregate(sizeblock~ngen, data= recom_block,mean)
means$sizeblock <- round(means$sizeblock, 2)

pdf(file = "data/Jackson_Lab_12_batches/boxplot_mean_recomb_block_size.pdf", height = 8, width = 10)
p1 <- ggplot(recom_block, aes(x=ngen, y=sizeblock, group = ngen, fill = ngen)) + 
  geom_boxplot(show.legend = F , outlier.size = 0.5, notchwidth = 3) +
  scale_x_discrete(drop=FALSE, breaks = c(21:23,NA,25,rep(NA,3),29:36)) +
  scale_fill_discrete_qualitative(palette = "warm")+
  geom_text(data = means, alpha = 0.85, aes(label = sizeblock, y = sizeblock + 0.15 )) + 
  ylab("Recombination Block Size (Mb)") +
  xlab("Generation") +
  labs(fill = "") +
  #ylim(c(0, 60)) +
  scale_y_continuous(breaks=c(0,5,10, 20, 40, 60), limits=c(0, 60)) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        text = element_text(size=16), 
        axis.title=element_text(size=16)) +
  guides(shape = guide_legend(override.aes = list(size = 12)))
p1
dev.off()

png 
  2 

p1

Past versions of “Average haplotype block size-1.png”

#block size distribution
for(g in unique(gm_after_qc$covar$ngen)){
  #plot for recom block size
  #png(paste0("data/Jackson_Lab_12_batches/DO_recom_block_size_G", g, ".png"))
  x <- pos_ind_gen[[g]][pos_ind_gen[[g]] != 0]
  # estimate the parameters
  fit1 <- fitdistr(x, "exponential") 
  
  # goodness of fit test
  ks.test(x, "pexp", fit1$estimate) # p-value > 0.05 -> distribution not refused
  
  # plot a graph
  hist(x, 
       freq = FALSE, 
       breaks = 200, 
       xlim = c(0, 5+quantile(x, 1)), 
       #ylim = c(0,0.3),
       xlab = "Recombination Block Size (Mb)", 
       main = paste0("Gen ", g))
  curve(dexp(x, rate = fit1$estimate), 
        from = 0, 
        to = 5+quantile(x, 1), 
        col = "red", 
        add = TRUE)
  #dev.off()
}

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-2.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-3.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-4.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-5.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-6.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-7.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-8.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-9.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-10.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-11.png”

Warning in ks.test(x, "pexp", fit1$estimate): ties should not be present for the
Kolmogorov-Smirnov test

Past versions of “Average haplotype block size-12.png”

Past versions of “Average haplotype block size-13.png”

#Average heterozygosity value

load("data/Jackson_Lab_12_batches/dat_het_ind_pr.RData")
dat_het_ind_pr$ngen <- factor(dat_het_ind_pr$ngen, levels = as.character(c(21:36)))

pdf(paste("data/Jackson_Lab_12_batches/DO_Heterozygosity_value_violin_genoprops.pdf"), width = 10, height =8)
p2 <- ggplot(dat_het_ind_pr, aes(x=ngen, y=het, group=ngen, fill=ngen)) +
  geom_violin(show.legend = FALSE) +
  geom_boxplot(show.legend = FALSE, width=0.35, color="black", alpha=0.6) +
  scale_x_discrete(drop=FALSE, breaks = c(21:23,NA,25,rep(NA,3),29:36)) +
  scale_fill_discrete_qualitative(palette = "warm")+
  ylab("Heterozygosity from genotype props") +
  xlab("Generation") +
  ylim(c(0.65, 1)) +
  geom_hline(yintercept=0.875, linetype="dashed", color = "red") +
  #scale_y_continuous(breaks=c(0.55, 0.65, 0.75, 0.85, 0.95, 1), limits=c(0.55, 1)) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        text = element_text(size=16),
        axis.title=element_text(size=16),
        legend.title=element_blank())
p2
dev.off()

png 
  2 

p2

Past versions of “Average heterozygosity value-1.png”

Session information

sessionInfo()

R version 4.0.0 (2020-04-24)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.4 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] colorspace_1.4-1  plotly_4.9.2.1    MASS_7.3-51.6     vcd_1.4-8        
 [5] ggplot2_3.3.2     lme4_1.1-23       Matrix_1.2-18     regress_1.3-21   
 [9] gap_1.2.2         abind_1.4-5       doParallel_1.0.15 iterators_1.0.12 
[13] foreach_1.5.0     data.table_1.12.8 dplyr_1.0.0       tidyr_1.1.0      
[17] table1_1.2        qtl2_0.22-8       workflowr_1.6.2  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6      lattice_0.20-41   zoo_1.8-8         rprojroot_1.3-2  
 [5] digest_0.6.25     lmtest_0.9-38     R6_2.4.1          backports_1.1.6  
 [9] RSQLite_2.2.0     evaluate_0.14     httr_1.4.1        pillar_1.4.4     
[13] rlang_0.4.6       lazyeval_0.2.2    minqa_1.2.4       whisker_0.4      
[17] nloptr_1.2.2.2    blob_1.2.1        rmarkdown_2.5     labeling_0.4.2   
[21] splines_4.0.0     statmod_1.4.34    stringr_1.4.0     htmlwidgets_1.5.1
[25] bit_1.1-15.2      munsell_0.5.0     compiler_4.0.0    httpuv_1.5.4     
[29] xfun_0.13         pkgconfig_2.0.3   htmltools_0.4.0   tidyselect_1.1.0 
[33] tibble_3.0.1      codetools_0.2-16  viridisLite_0.3.0 crayon_1.3.4     
[37] withr_2.2.0       later_1.0.0       jsonlite_1.6.1    nlme_3.1-147     
[41] gtable_0.3.0      lifecycle_0.2.0   DBI_1.1.0         git2r_0.27.1     
[45] magrittr_1.5      scales_1.1.1      stringi_1.4.6     farver_2.0.3     
[49] fs_1.4.1          promises_1.1.0    ellipsis_0.3.0    generics_0.0.2   
[53] vctrs_0.3.1       boot_1.3-25       Formula_1.2-4     tools_4.0.0      
[57] bit64_0.9-7       glue_1.4.0        purrr_0.3.4       yaml_2.2.1       
[61] memoise_1.1.0     knitr_1.28