inst/doc/NIPTeR.R

In NIPTeR: Fast and Accurate Trisomy Prediction in Non-Invasive Prenatal Testing

## ---- echo = FALSE-------------------------------------------------------
knitr::opts_chunk$set(
  eval = FALSE
  
)

## ---- eval = TRUE, echo = FALSE------------------------------------------
example_bin_counts <- structure(c(308, 243, 0, 624, 448, 0, 247, 163, 303, 263, 189, 
637, 249, 174, 266, 216, 270, 262, 287, 257, 274, 213, 234, 244, 
157, 269, 245, 240, 262, 142, 217, 209, 148, 225, 264, 229, 239, 
231, 142, 243, 99, 225, 251, 241, 205, 324, 88, 255, 146, 285, 
258, 251, 241, 355, 264, 259, 236, 259, 242, 216, 244, 383, 230, 
261), .Dim = c(8L, 8L), .Dimnames = list(c(" 1", "2", "3", "4", 
"5", "6", "7", "8"), c("1", "2", "3", "4", "5", "6", "7", "8"
)))

## ---- eval = TRUE, echo = FALSE------------------------------------------
example_regression_stats <- structure(c("-0.171248675025602", "0.00615122686158584", "Practical_CV", 
"0.583158688871931", "11  2  5  12", "1.00050469190991", "0.00335713985926308", 
"-0.799303637284865", "0.00599131304613039", "Practical_CV", 
"0.503865450973465", "10  19  17  16", "1.00011356353248", "0.00394053718934269"
), .Dim = c(7L, 2L), .Dimnames = list(c("Z_score_sample", "CV", 
"cv_types", "P_value_shapiro", "Predictor_chromosomes", "Mean_test_set", 
"CV_train_set"), c("Prediction_set_1", "Prediction_set_2")))

## ---- eval = TRUE, echo = FALSE------------------------------------------
example_bin_counts_sex <- structure(c(0, 0, 160, 0, 6, 0, 266, 0, 127, 0, 85, 0, 182, 0, 
235, 0), .Dim = c(2L, 8L), .Dimnames = list(c("X", "Y"), c("1", 
"2", "3", "4", "5", "6", "7", "8")))

## ---- eval = TRUE, echo = FALSE------------------------------------------
example_match_report <- structure(c(6.76557504762928e-07, 2.15862500903867e-06, 3.10100948006486e-06, 
3.29198583437393e-06, 3.58809416090935e-06, 3.83885934221371e-06, 
4.43239735020721e-06, 7.8916102881966e-06, 9.73718696931918e-06, 
1.12605046722957e-05), .Dim = c(10L, 1L), .Dimnames = list(c("Sample 6", 
"Sample 8", "Sample 1", "Sample 10", "Sample 9", "Sample 3", "Sample 5", 
"Sample 2", "Sample 7", "Sample 4"), "Sum_of_squares"))

## ---- eval=TRUE, echo = FALSE, results = "asis"--------------------------
pander::pandoc.table(example_bin_counts, caption = "Table 1: Example data bin counts. The rows represent chromosomes and the columns bins", include.rownames = T)

## ---- eval=TRUE, echo = FALSE, results = "asis"--------------------------
pander::pandoc.table(example_bin_counts_sex, caption = "Table 2: Example data bin counts sex chromosomes. The rows represent chromosomes and the columns bins", include.rownames = T)

## ------------------------------------------------------------------------
#  bin_bam_sample(bam_filepath, do_sort = F, separate_strands = F, custom_name = NULL)

## ------------------------------------------------------------------------
#  gc_correct(nipt_object, method, include_XY = F, span = 0.75)

## ------------------------------------------------------------------------
#  #Correct NIPTSample object using LOESS method
#  loess_corrected_sample <- gc_correct(nipt_object = sample_of_interest, method = "LOESS",
#                                       include_XY = F, span = 0.75)
#  #Correct NIPTControlGroup object using bin method
#  gc_bin_corrected_control_group <- gc_correct(nipt_object = control_group, method = "bin",
#                                               include_XY = T)
#  

## ------------------------------------------------------------------------
#  matchcontrolgroup(nipt_sample, nipt_control_group, mode, n_of_samples,
#    include_chromosomes = NULL, exclude_chromosomes = NULL)

## ---- eval=TRUE, echo = FALSE, results = "asis"--------------------------
pander::pandoc.table(example_match_report , caption = "Table 3: Example match_control_group mode 'subset'", include.rownames = T)

## ------------------------------------------------------------------------
#  #Run matchcontrolgroup in mode "report"
#  scores_control_group <- matchcontrolgroup(nipt_sample = gc_LOESS_corrected_sample,
#                                            nipt_control_group = gc_LOESS_corrected_control_group,
#                                            mode = "report", include_chromosomes = c(13,18))
#  #Run matchcontrolgroup in mode "subset" and select 50 best matching samples
#  subset_loess_corrected_control_group <- matchcontrolgroup(nipt_sample = gc_LOESS_corrected_sample,
#                                                            nipt_control_group = loess_corrected_control_group,
#                                                            mode = "subset", n_of_samples = 50)
#  

## ------------------------------------------------------------------------
#  chi_correct(nipt_sample, nipt_control_group, chi_cutoff = 3.5, include_XY = F)

## ------------------------------------------------------------------------
#  #Apply chi-squared based variation reduction method
#  chi_corrected_data <- chicorrect(nipt_sample = gc_LOESS_corrected_sample,
#                                   nipt_control_group = subset_loess_corrected_control_group)
#  #Extract sample and control group
#  loess_chi_corrected_sample <- chi_corrected_data$sample
#  subset_loess_chi_corrected_control_group <- chi_corrected_data$control_group

## ------------------------------------------------------------------------
#  calculate_z_score(nipt_sample, nipt_control_group, chromo_focus)

## ------------------------------------------------------------------------
#  #Calculate Z score for chromosome 13
#  z_score_result_13 <- calculate_z_score(nipt_sample = loess_chi_corrected_sample,
#                                         nipt_control_group = subset_loess_chi_corrected_control_group,
#                                         chromo_focus = 13)

## ------------------------------------------------------------------------
#  perform_regression(nipt_sample, nipt_control_group, chromo_focus,
#    n_models = 4, n_predictors = 4, exclude_chromosomes = NULL,
#    include_chromosomes = NULL, use_test_train_set = T,
#    size_of_train_set = 0.6, overdispersion_rate = 1.15,
#    force_practical_vc = F)

## ---- eval=TRUE, echo = FALSE, results = "asis"--------------------------
pander::pandoc.table(example_regression_stats, caption = "Table 2: Example regression results and statistics", include.rownames = T)

## ------------------------------------------------------------------------
#  prepare_ncv(nipt_control_group, chr_focus, max_elements,
#    exclude_chromosomes = NULL, include_chromosomes = NULL,
#    use_test_train_set = T, size_of_train_set = 0.6)

## ------------------------------------------------------------------------
#  calculate_ncv_score(nipt_sample, ncv_template)

## ------------------------------------------------------------------------
#  as_control_group(nipt_samples)

## ------------------------------------------------------------------------
#  add_samples_controlgroup(nipt_control_group, samples_to_add)

## ------------------------------------------------------------------------
#  remove_duplicates_controlgroup(nipt_control_group)

## ------------------------------------------------------------------------
#  remove_sample_controlgroup(samplename, nipt_control_group)

## ------------------------------------------------------------------------
#  diagnose_control_group(nipt_control_group)

## ------------------------------------------------------------------------
#  
#  list_NIPT_samples <- lapply(X = bam_filepaths, bin_bam_sample, do_sort = FALSE,
#                               separate_strands = T, custom_name = NULL)
#  

## ------------------------------------------------------------------------
#  control_group  <- as_control_group(nipt_samples = lapply(X = bam_filepaths, bin_bam_sample,
#                                                           do_sort = F, separate_strands = T))

## ------------------------------------------------------------------------
#  control_group  <- as_control_group(nipt_samples = mapply(FUN = bin_bam_sample, bam_filepaths,
#                                                           custom_name = names, SIMPLIFY = FALSE))

## ------------------------------------------------------------------------
#  library(NIPTeR)
#  #Retrieve filenames
#  bam_filepaths <- list.files(path = "/Path/to/bamfiles/", pattern = ".bam", full.names = T)
#  #Load files and convert to control group
#  control_group  <- as_control_group(nipt_samples = lapply(X = bam_filepaths, bin_bam_sample,
#                                                           do_sort = F, separate_strands = FALSE))
#  #Save control group for later
#  saveRDS(object = control_group, file = "/Path/to/directory/control_group.rds")

## ------------------------------------------------------------------------
#  library(NIPTeR)
#  #Gather all bam filepaths in a vector. Corresponds to 1a in figure
#  bam_filepaths <- list.files(path = "/Path/to/bamfiles/", pattern = ".bam", full.names = T)
#  #Load all bam files using lapply and feed the results to function as_control_group,
#  #converting the NIPTSamples to a NIPTControlGroup object. Corresponds to 1b in figure
#  control_group  <- as_control_group(nipt_samples = lapply(X = bam_filepaths, bin_bam_sample,
#                                                           do_sort = F, separate_strands = F))
#  #apply a gc LOESS correction to all samples. Since this can take up to 30 seconds
#  #sample, doing this once for a control group and store the results can be rewarding
#  #in terms of analysis time. Corresponds to 2a in figure
#  gc_control_group <- gc_correct(nipt_object = control_group, method = "LOESS")
#  #Retrieve control group diagnostics. Corresponds with 3a in figure
#  control_group_diagnostics <- diagnose_control_group(nipt_control_group = control_group)
#  #Retrieve samplenames with an abberant Z score for any chromosome and remove these samples
#  #from the control group. Corresponds with 3b in figure
#  abberant_sample_names <- unique(control_group_diagnostics$abberant_scores$Sample_name)
#  for (i in 1:length(abberant_sample_names)){
#    control_group <- remove_sample_controlgroup(samplename = abberant_sample_names[i],
#                                                nipt_control_group = control_group)
#  }
#  #Save the gc_corrected control groups to disk. Corresponds to 4a in figure
#  saveRDS(object = control_group, file = "/path/to/controlgroup.rds")

## ------------------------------------------------------------------------
#  library(NIPTeR)
#  #Load sample. Corresponds with 1a in figure
#  sample_of_interest <- bin_bam_sample(bam_filepath = "/Path/to/bamfile.bam", separate_strands = T)
#  #Load control group. Corresponds with 1b in figure
#  control_group <- readRDS("/Path/to/control_group.rds")
#  #Perform a chi-square based variation reduction on new trimmed control group and sample and
#  #extract data. Corresponds with 2a in figure
#  chi_data <- chi_correct(nipt_sample = sample_of_interest, nipt_control_group = control_group)
#  sample_of_interest <- chi_data$sample
#  control_group <- chi_data$control_group
#  #Perform regression for chromosome 21 with default settings, so:
#  #Create 4 models with 4 predictors each
#  #All chromosomes are potential predictors except the potential trisomic chromosomes 13, 18 and 21
#  #Use a test and train set where the size of the train set is 60% of the control group
#  #Assume at least 15% overdispersion in the data
#  #Dont force practical CV, so if the CV is below 1.15 * standard error of the mean use this as VC
#  #Corresponds with 3a in figure
#  regression_score_21 <- perform_regression(nipt_sample = sample_of_interest,
#                                            nipt_control_group = control_group, chromo_focus = 21)
#  
#  
#  ###       Gather relevant data from the objects on the workspace     #######

## ------------------------------------------------------------------------
#  library(NIPTeR)
#  #Load sample. Corresponds with 1a in figure
#  sample_of_interest <- bin_bam_sample(bam_filepath = "/Path/to/bamfile.bam")
#  #Load control group. Corresponds with 1b in figure
#  control_group <- readRDS("/Path/to/control_group.rds")
#  #Peform a GC correction type bin for the sample and the controlgroup
#  #Corresponds with 2a in figure
#  sample_of_interest <- gc_correct(nipt_object = sample_of_interest, method = "bin")
#  control_group <- gc_correct(nipt_object = control_group, method = "bin")
#  #Trim control group by only selecting 80% of best matching control samples
#  #Corresponds with 2b in figure
#  control_group <- match_control_group(nipt_sample = sample_of_interest, nipt_control_group = control_group,
#                                       mode = "subset", n_of_samples = round(length(control_group$samples) *.8,
#                                                                             digits = 0))
#  #Perform a chi-square based variation reduction on new trimmed control group and sample and
#  #extract data. Corresponds with 2c in figure
#  chi_data <- chi_correct(nipt_sample = sample_of_interest, nipt_control_group = control_group)
#  sample_of_interest <- chi_data$sample
#  control_group <- chi_data$control_group
#  #Retrieve control group diagnostics. Corresponds with 3a in figure
#  control_group_diagnostics <- diagnose_control_group(nipt_control_group = control_group)
#  #Retrieve samplenames with an abberant Z score for any chromosome and remove these samples
#  #from the control group. Corresponds with 3b in figure
#  abberant_sample_names <- unique(control_group_diagnostics$abberant_scores$Sample_name)
#  for (i in 1:length(abberant_sample_names)){
#    control_group <- remove_sample_controlgroup(samplename = abberant_sample_names[i],
#                                                nipt_control_group = control_group)
#  }
#  #Calculate Z score from chromosomes 13, 18 and 21. Corresponds with 4a in figure
#  z_score_13 <- calculate_z_score(nipt_sample = sample_of_interest,
#                                  nipt_control_group = control_group, chromo_focus = 13)
#  z_score_18 <- calculate_z_score(nipt_sample = sample_of_interest,
#                                  nipt_control_group = control_group, chromo_focus = 18)
#  z_score_21 <- calculate_z_score(nipt_sample = sample_of_interest,
#                                  nipt_control_group = control_group, chromo_focus = 21)
#  #Perform regression for all potential trisomic chromosomes with default settings, so:
#  #Create 4 models for every potential trisomic chromosome with 4 predictors each
#  #All chromosomes are potential predictors except the potential trisomic chromosomes 13, 18 and 21
#  #Use a test and train set where the size of the train set is 60% of the control group
#  #Assume at least 15% overdispersion in the data
#  #Dont force practical CV, so if the CV is below 1.15 * standard error of the mean use this as VC
#  #Corresponds with 4c in figure
#  regression_score_13 <- perform_regression(nipt_sample = sample_of_interest,
#                                            nipt_control_group = control_group, chromo_focus = 13)
#  regression_score_18 <- perform_regression(nipt_sample = sample_of_interest,
#                                            nipt_control_group = control_group, chromo_focus = 18)
#  regression_score_21 <- perform_regression(nipt_sample = sample_of_interest,
#                                            nipt_control_group = control_group, chromo_focus = 21)
#  #Get NCVTemplates for all potential trisomic chromosomes with max 9 denominators and default settings, so:
#  #All autosomals chromosomes are potential predictors, except the potential trisomic chromosomes 13, 18 and 21
#  #Use a test and train set where the size of the train set is 60% of the control group
#  #Corresponds with 4c in figure
#  new_ncv_template_13 <- prepare_ncv(nipt_control_group = control_group, chr_focus = 13, max_elements = 9)
#  new_ncv_template_18 <- prepare_ncv(nipt_control_group = control_group, chr_focus = 18, max_elements = 9)
#  new_ncv_template_21 <- prepare_ncv(nipt_control_group = control_group, chr_focus = 21, max_elements = 9)
#  #Use the NCVTemplates to get NCV scores for the sample of interest
#  #Corresponds with 4d in figure
#  ncv_score_13 <- calculate_ncv_score(nipt_sample = sample_of_interest, ncv_template = new_ncv_template_13)
#  ncv_score_18 <- calculate_ncv_score(nipt_sample = sample_of_interest, ncv_template = new_ncv_template_18)
#  ncv_score_21 <- calculate_ncv_score(nipt_sample = sample_of_interest, ncv_template = new_ncv_template_21)
#  
#  
#  ###       Gather relevant data from the objects on the workspace     #######
#  
#  
#