## ----dependencies, warning=FALSE, message=FALSE------------------------------- library(bsseq) library(tidyverse) #Additional software modules needed for preprocessing: #Dorado https://github.com/nanoporetech/dorado #SAMtools https://github.com/samtools #Minimap2 https://github.com/lh3/minimap2 #Modkit https://github.com/nanoporetech/modkit ## ----read all, message=TRUE--------------------------------------------------- files <- list.files(system.file("extdata/longread/silversides/all_cpgs", package = "bsseq"), full.names = TRUE) mc_all<-read.bedMethyl(files, strandCollapse = T, output = "MethylCounts") mc_all ## ----read cpg, message=FALSE-------------------------------------------------- files <- list.files(system.file("extdata/longread/silversides/ref_cpgs", package = "bsseq"), full.names = TRUE) mc_cpg<-read.bedMethyl(files=files, strandCollapse = T, output = "MethylCounts") mc_cpg ## ----cov filtering all single sample------------------------------------------ #get the first sample mc_all_1 <- mc_all[,1] #get the indices of loci with coverage >= 5 loci.idx <- which(getMethylCounts(mc_all_1, type="Cov")>= 5) #filter the object to retain these loci mc_all_1_filtered <- mc_all_1[loci.idx,] mc_all_1_filtered ## ----cov filtering cpg single sample------------------------------------------ mc_cpg_1 <- mc_cpg[,1] loci.idx <- which(getMethylCounts(mc_cpg_1, type="Cov")>= 5) mc_cpg_1_filtered<-mc_cpg_1[loci.idx,] mc_cpg_1_filtered ## ----cov filtering all and cpg single sample, echo=FALSE---------------------- results <- tibble() for (i in 1:30) { sites_all_count <- length(which(DelayedMatrixStats::rowSums2(getMethylCounts( mc_all[,1], type = "Cov") >= i) == 1)) sites_cpg_count <- length(which(DelayedMatrixStats::rowSums2(getMethylCounts( mc_cpg[,1], type = "Cov") >= i) == 1)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count, ref_CpG = sites_cpg_count)) } #make results long format results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = threshold, y= sites, color = Type)) + geom_point() + xlab("Coverage threshold") + ylab("Number of loci") + ggtitle("Coverage filtering") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----.99 all filtering all---------------------------------------------------- #get the first sample mc_all_1 <- mc_all[,1] #get the indices of loci with scaled likelihood above 0.99 of being a "allCpG" loci.idx <- getCpGs(mc_all_1, type = "allCpG", threshold = 0.99) #filter the object to retain these loci mc_all_1_filtered <- mc_all_1[loci.idx,] mc_all_1_filtered ## ----.99 all filtering cpg---------------------------------------------------- mc_cpg_1 <- mc_cpg[,1] loci.idx <- getCpGs(mc_cpg_1, type = "allCpG", threshold = 0.99) mc_cpg_1_filtered <- mc_cpg_1[loci.idx,] mc_cpg_1_filtered ## ----all likelihood threshold test, echo=FALSE-------------------------------- results <- tibble() for (i in seq(0.95, 1, 0.0005)) { sites_all_count <- length( getCpGs(mc_all[,1], type = "allCpG", threshold = i)) sites_cpg_count <- length( getCpGs(mc_cpg[,1], type = "allCpG", threshold = i)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count, ref_CpG = sites_cpg_count)) } #change the code to run from 0.05 to 1 with 0.05 increments results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = threshold, y= sites, color = Type)) + geom_point() + xlab("Likelihood threshold") + ylab("Number of loci") + ggtitle("Homozygous and heterozygous CpG loci") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----.99 hom filtering all---------------------------------------------------- #get the first sample mc_all_1 <- mc_all[,1] #get the indices of loci with scaled likelihood above 0.99 of being homozygous loci.idx <- getCpGs(mc_all_1, type = "homozygous", threshold = 0.99) #filter the object to retain these loci mc_all_1_filtered <- mc_all_1[loci.idx,] mc_all_1_filtered ## ----.99 hom filtering cpg---------------------------------------------------- mc_cpg_1 <- mc_cpg[,1] loci.idx <- getCpGs(mc_cpg_1, type = "homozygous", threshold = 0.99) mc_cpg_1_filtered <- mc_cpg_1[loci.idx,] mc_cpg_1_filtered ## ----hom likelihood threshold test, echo=FALSE-------------------------------- results <- tibble() for (i in seq(0.95, 1, 0.0005)) { sites_all_count_hom <- length( getCpGs(mc_all[,1], type = "homozygous", threshold = i)) sites_cpg_count_hom <- length( getCpGs(mc_cpg[,1], type = "homozygous", threshold = i)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count_hom, ref_CpG = sites_cpg_count_hom)) } #change the code to run from 0.05 to 1 with 0.01 increments results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = threshold, y= sites, color = Type)) + geom_point() + xlab("Likelihood threshold") + ylab("Number of loci") + ggtitle("Homozygous CpG loci") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----.99 het filtering all---------------------------------------------------- #get the first sample mc_all_1 <- mc_all[,1] #get the indices of loci with scaled likelihood above 0.99 of being heterozygous loci.idx <- getCpGs(mc_all_1, type = "heterozygous", threshold = 0.99) #filter the object to retain these loci mc_all_1_filtered <- mc_all_1[loci.idx,] mc_all_1_filtered ## ----.99 het filtering cpg---------------------------------------------------- mc_cpg_1 <- mc_cpg[,1] loci.idx <- getCpGs(mc_cpg_1, type = "heterozygous", threshold = 0.99) mc_cpg_1_filtered <- mc_cpg_1[loci.idx,] mc_cpg_1_filtered ## ----het likelihood threshold test, echo=FALSE-------------------------------- results <- tibble() for (i in seq(0.95, 1, 0.0005)) { sites_all_count_het <- length( getCpGs(mc_all[,1], type = "heterozygous", threshold = i)) sites_cpg_count_het <- length( getCpGs(mc_cpg[,1], type = "heterozygous", threshold = i)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count_het, ref_CpG = sites_cpg_count_het)) } #change the code to run from 0.05 to 1 with 0.01 increments results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = threshold, y= sites, color = Type)) + geom_point() + xlab("Likelihood threshold") + ylab("Number of loci") + ggtitle("Heterozygous CpG loci") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----multi coverage threshold, echo=FALSE------------------------------------- results <- tibble() for (i in 1:30) { sites_all_count <- length(which( DelayedMatrixStats::rowSums2(getMethylCounts(mc_all, type = "Cov") >= i) == 9)) sites_cpg_count <- length(which( DelayedMatrixStats::rowSums2(getMethylCounts(mc_cpg, type = "Cov") >= i) == 9)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count, ref_CpG = sites_cpg_count)) } #make results long format results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = as.factor(threshold), y= sites, color = Type)) + geom_point() + xlab("Coverage threshold") + ylab("Number of loci") + ggtitle("Coverage filtering all samples") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----multi samples threshold, echo=FALSE-------------------------------------- results <- tibble() for (i in 1:9) { sites_all_count <- length(which(DelayedMatrixStats::rowSums2(getMethylCounts( mc_all, type = "Cov") >= 5) >= i)) sites_cpg_count <- length(which(DelayedMatrixStats::rowSums2( getMethylCounts(mc_cpg, type = "Cov") >= 5) >= i)) results <- bind_rows(results, tibble(threshold = i, all_CpG = sites_all_count, ref_CpG = sites_cpg_count)) } #make results long format results_long <- results %>% pivot_longer(cols = c(all_CpG, ref_CpG), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = as.factor(threshold), y= sites, color = Type)) + geom_point() + xlab("Samples passing 5X coverage threshold") + ylab("Number of loci") + ggtitle("Coverage filtering") + scale_color_manual(values = c("all_CpG" = "blue", "ref_CpG" = "red")) + theme_bw(base_size = 14) ## ----multi coverage and samples filtering------------------------------------- loci.idx <- which( DelayedMatrixStats::rowSums2(getMethylCounts(mc_all, type="Cov")>= 5) >= 6) mc_coverage_filtered <- mc_all[loci.idx,] mc_coverage_filtered ## ----getCpGMatrix------------------------------------------------------------- G_all <- getCpGMatrix(mc_all) head(G_all) ## ----getMaxLikelihoodMatrix--------------------------------------------------- Q_all <- getMaxLikelihoodMatrix(mc_all) head(round(Q_all,3)) ## ----getCpGMatrix allCpG------------------------------------------------------ G_all <- getCpGMatrix(mc_all, allCpG = TRUE) head(G_all) ## ----getMaxLikelihoodMatrix allCpG-------------------------------------------- Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE) head(round(Q_all,3)) ## ----multi likelihood threshold test all, echo=FALSE-------------------------- G_all <- getCpGMatrix(mc_all, allCpG = TRUE) G_cpg <- getCpGMatrix(mc_cpg, allCpG = TRUE) Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE) Q_cpg <- getMaxLikelihoodMatrix(mc_cpg, allCpG = TRUE) results_Q <- tibble() for (i in seq(0.95, 1, 0.0005)) { sites_count_all <- length(which( DelayedMatrixStats::rowSums2(G_all == 0) == 9 & DelayedMatrixStats::rowSums2(Q_all >= i) == 9)) sites_count_cpg <- length(which( DelayedMatrixStats::rowSums2(G_cpg == 0) == 9 & DelayedMatrixStats::rowSums2(Q_cpg >= i) == 9)) results_Q <- bind_rows(results_Q, tibble(threshold = i, all_CpG_hom = sites_count_all, ref_CpG_hom = sites_count_cpg)) } #make results long format results_long <- results_Q %>% pivot_longer(cols = c(all_CpG_hom, ref_CpG_hom), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = threshold, y= sites, color = Type)) + geom_point() + xlab("Likelihood threshold") + ylab("Number of CpG loci") + ggtitle("Likelihood filtering") + scale_color_manual(values = c("all_CpG_hom" = "blue", "ref_CpG_hom" = "red")) + theme_bw(base_size = 14) ## ----multi samples threshold test all, echo=FALSE----------------------------- G_all <- getCpGMatrix(mc_all, allCpG = TRUE) G_cpg <- getCpGMatrix(mc_cpg, allCpG = TRUE) Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE) Q_cpg <- getMaxLikelihoodMatrix(mc_cpg, allCpG = TRUE) results <- tibble() for (i in 1:9) { sites_count_all <- length(which( DelayedMatrixStats::rowSums2(G_all == 0 & Q_all >= 0.99 ) >= i)) sites_count_cpg <- length(which( DelayedMatrixStats::rowSums2(G_cpg == 0 & Q_cpg >= 0.99 ) >= i)) results <- bind_rows(results, tibble(threshold = i, all_CpG_hom = sites_count_all, ref_CpG_hom = sites_count_cpg)) } #make results long format results_long <- results %>% pivot_longer(cols = c(all_CpG_hom, ref_CpG_hom), names_to = "Type", values_to = "sites") #plot the results with ggplot(results_long, aes(x = as.factor(threshold), y= sites, color = Type)) + geom_point() + xlab("Samples passing .99 all filtering") + ylab("Number of CpG loci") + ggtitle("Samples passing likelihood filtering") + scale_color_manual(values = c("all_CpG_hom" = "blue", "ref_CpG_hom" = "red"))+ theme_bw(base_size = 14) ## ----multi likelihood and samples filtering----------------------------------- G_all <- getCpGMatrix(mc_all, allCpG = TRUE) Q_all <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE) loci.idx <- which( DelayedMatrixStats::rowSums2(Q_all >= .99 & G_all==0) >= 6) mc_filtered <- mc_all[loci.idx,] mc_filtered ## ----------------------------------------------------------------------------- #hdf5_path <- # "~/Documents/bsseq/inst/extdata/longread/silversides/mc_filtered.hdf5" #HDF5Array::saveHDF5SummarizedExperiment(mc_filtered, hdf5_path) # Load pre-saved HDF5 file files <- system.file("extdata/longread/silversides/mc_filtered.hdf5", package = "bsseq") mc_filtered <- HDF5Array::loadHDF5SummarizedExperiment(files) mc_filtered ## ----------------------------------------------------------------------------- # Filter and process a sample from the MethylCounts object mc_sample <- mc_filtered[, 1] mc_sample_filtered <- mc_sample[ getCpGs(mc_sample, type = "homozygous", threshold = 0.99)] # Display the filtered MethylCounts object mc_sample_filtered ## ----------------------------------------------------------------------------- # Convert to BSseq with both 5mC and 5hmC data bs <- BSseq(mc = mc_sample_filtered, mods = "5mC+5hmC") bs # Calculate the mean methylation using raw data mean(getMeth(bs, type = "raw")) ## ----------------------------------------------------------------------------- # Convert to BSseq with only 5mC data bs_M <- BSseq(mc = mc_sample_filtered, mods = "5mC") bs_M # Calculate the mean methylation for 5mC mean(getMeth(bs_M, type = "raw")) ## ----------------------------------------------------------------------------- # Convert to BSseq with only 5hmC data bs_H <- BSseq(mc = mc_sample_filtered, mods = "5hmC") bs_H # Calculate the mean methylation for 5hmC mean(getMeth(bs_H, type = "raw")) ## ----multi coverage samples and likelihood filtering, message=FALSE----------- files <- list.files(system.file("extdata/longread/silversides/all_cpgs", package = "bsseq"), full.names = TRUE) mc_all<-read.bedMethyl(files, strandCollapse = T, output = "MethylCounts") G <- getCpGMatrix(mc_all, allCpG = TRUE) Q <- getMaxLikelihoodMatrix(mc_all, allCpG = TRUE) loci.idx <- which( DelayedMatrixStats::rowSums2(Q >= .99 & G==0) >= 6 & DelayedMatrixStats::rowSums2(getMethylCounts(mc_all, type="Cov")>= 5) >= 6) mc_filtered <- mc_all[loci.idx,] mc_filtered ## ----coverage samples and likelihood filtering-------------------------------- mc_sample<- mc_filtered[,1] mc_sample_filtered<-mc_sample[ getCpGs(mc_sample, type = "homozygous", threshold = 0.99)] mc_sample_filtered ## ----------------------------------------------------------------------------- bs_sample_filtered<-BSseq(mc = mc_sample_filtered) bs_sample_filtered ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()