R Package MADSEQ

In MADSEQ: Mosaic Aneuploidy Detection and Quantification using Massive Parallel Sequencing Data

The MADSEQ package is a group of hierarchical Bayesian model for the detection and quantification of potential mosaic aneuploidy in sample using massive parallel sequencing data.

The MADSEQ package takes two pieces of information for the detection and quantification of mosaic aneuploidy:

1. The distribution of the alternative allele frequencies (AAF) of the sites that are genotyped as heterozygous.
2. The average sequencing coverage for regions. (for targeted sequencing it's each targeted region; for whole genome sequencing, bin the genome into bins. And because sequencing coverage are usually biased by GC content, normalization is necessary, the normalization function is provided in the package).

MADSEQ works on the whole chromosome resolution. It applies all of the five models (normal, monosomy, mitotic trisomy, meiotic trisomy, loss of heterozygosity) to fit the distribution of the AAF of all the heterozygous sites, and fit the distribution of the coverage from that chromosome. After fitting the same data using all models, it does model comparison using BIC (Bayesian Information Criteria) to select the best model. The model selected tells us whether the chromosome is aneuploid or not, and also the type of mosaic aneuploidy. Then, from the posterior distribution of the best model, we could get the estimation of the fraction of aneuploidy cells.

Data You Need

• A bam file for your sample: sorted, indexed.
• A vcf.gz file for your sample:
  • The vcf should be bgzip compressed (bgzip is a part of tabix package, which can be download from here).
  • The FORMAT section of the vcf file should contain depth (DP) and allelic depth (AD) field.
• A bed file containing targeted regions: if it's targeted sequencing data, then the bed file contatining targeted regions (tiled regions is recommended) . If it's whole genome sequencing, bin the whole genome into equal size bins (i.e. 500kb).

Note: Currently our package only supports one bam and one vcf file per sample. If you have more than one sample, please prepare multiple bam and vcf files for each of them.

Example Data

There are two sets of example data come with the package:

1. A bam file named aneuploidy.bam and a vcf file named aneuploidy.vcf.gz for a sample containing trisomy in chromosome 18.
2. A bam file named normal.bam and a vcf file names normal.vcf.gz for a normal sample.
3. A bed file named target.bed containing the targeted regions.

knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

To access the data use

system.file("extdata","aneuploidy.bam",package="MADSEQ")
system.file("extdata","aneuploidy.vcf.gz",package="MADSEQ")

Note:This is just a set of example data, only contains a very little region of the genome.

Use The Package

We will start with the bam file, vcf file and bed file in the example data to show you each step for the analysis.

Prepare coverage information

Started with bam file and bed file, you can use prepareCoverageGC function to get the coverage and GC information for each targeted regions.

## load the package
suppressMessages(library("MADSEQ"))

## get path to the location of example data
aneuploidy_bam = system.file("extdata","aneuploidy.bam",package="MADSEQ")
normal_bam = system.file("extdata","normal.bam",package="MADSEQ")
target = system.file("extdata","target.bed",package="MADSEQ")

## Note: for your own data, just specify the path to the location 
## of your file using character.

## prepare coverage and GC content for each targeted region
# aneuploidy sample
aneuploidy_cov = prepareCoverageGC(target_bed=target, 
                                    bam=aneuploidy_bam, 
                                    "hg19")

# normal sample
normal_cov = prepareCoverageGC(target_bed=target, 
                                bam=normal_bam, 
                                "hg19")

## view the first two rows of prepared coverage data (A GRanges Object)
aneuploidy_cov[1:2]

normal_cov[1:2]

Normalize coverage information

The normalization function takes prepared coverage GRanges object from prepareCoverageGC function, normalize the coverage and calculate the expected coverage for the sample. If there is only one sample, the function will correct the coverage by GC content, and take the average coverage for the whole genome as expected coverage. If there are more than one samples given, the function will first quantile normalize coverage across samples, then correct the coverage by GC for each sample. If control sample is not specified, the expected coverage is the median coverage across all samples, if a normal control is specified, the average coverage for control sample is taken as expected coverage for further analysis.

Note:

1. If you have more than one samples, please make sure they are targeted by the same regions.
2. If you have more than one samples, please separate female samples and male samples, which means normalize them separately; because sex chromosomes can bias the normalization.
3. One sample is good to go, however multiple samples will help you get a better normalization.

If you only have one sample

If you choose to write the output to file (recommended)

## normalize coverage data
## set plot=FALSE here because similar plot will show in the following example
normalizeCoverage(aneuploidy_cov,writeToFile=TRUE, destination=".",plot=FALSE)

If you don't want to write output to file

## normalize coverage data
aneuploidy_normed = normalizeCoverage(aneuploidy_cov,writeToFile=FALSE,
                                      plot=FALSE)

## a GRangesList object will be produced by the function, look at it by
names(aneuploidy_normed)
aneuploidy_normed[["aneuploidy_cov"]]

If you have more than one samples, without a normal control

If you choose to write the output to file (recommended)

## normalize coverage data
normalizeCoverage(aneuploidy_cov, normal_cov,
                  writeToFile =TRUE, destination = ".", plot=FALSE)

If you don't want to write output to file

## normalize coverage data
normed_without_control = normalizeCoverage(aneuploidy_cov, normal_cov, 
                                           writeToFile=FALSE, plot=TRUE)

## a GRangesList object will be produced by the function
length(normed_without_control)
names(normed_without_control)

## subsetting
normed_without_control[["aneuploidy_cov"]]
normed_without_control[["normal_cov"]]

If you have more than one samples, with a normal control

If you choose to write the output to file (recommended)

## normalize coverage data, normal_cov is the control sample
normalizeCoverage(aneuploidy_cov, control=normal_cov,
                  writeToFile=TRUE, destination = ".",plot=FALSE)

If you don't want to write output to file

normed_with_control = normalizeCoverage(aneuploidy_cov, control=normal_cov, 
                                        writeToFile =FALSE, plot=FALSE)

## a GRangesList object will be produced by the function
length(normed_without_control)
names(normed_with_control)

Prepare heterozygous sites

Having vcf.gz file and target bed file ready, use prepareHetero function to process the heterozygous sites.

## specify the path to vcf.gz file
aneuploidy_vcf = system.file("extdata","aneuploidy.vcf.gz",package="MADSEQ")

## target bed file specified before

## If you choose to write the output to file (recommended)
prepareHetero(aneuploidy_vcf, target, genome="hg19", 
              writeToFile=TRUE, destination=".",plot = FALSE)

## If you don't want to write output to file
aneuploidy_hetero = prepareHetero(aneuploidy_vcf, target,
                                  genome="hg19", writeToFile=FALSE,plot = FALSE)

Run MadSeq model to detect potential mosaic aneuploidy

The function runMadSeq will run the models and select the best model for the input data.

Note:

1. Among three normalized coverage sets listed above (one sample, two samples without control, two samples with a control), we will use the two samplew with a control case for the following analysis.
2. Because these models are based on MCMC sampling, the running process can be very long. Running different chromosomes parallel in background or High Performance Computer Cluster is highly recommended.

If the processed data have been written into files

## specify the path to processed files
aneuploidy_hetero = "./aneuploidy.vcf.gz_filtered_heterozygous.txt"
aneuploidy_normed_cov = "./aneuploidy_cov_normed_depth.txt"

## run the model
aneuploidy_chr18 = runMadSeq(hetero=aneuploidy_hetero, 
                             coverage=aneuploidy_normed_cov, 
                             target_chr="chr18",
                             nChain=1, nStep=1000, thinSteps=1,
                             adapt=100,burnin=200)

## An MadSeq object will be returned
aneuploidy_chr18

Note: In order to save time, we only run 1 chain with a much less steps compared with default settings. For real cases, the default settings are recommended.

If the processed data have NOT been written into files

## subset normalized coverage for aneuploidy sample from the GRangesList 
## returned by normalizeCoverage function
aneuploidy_normed_cov = normed_with_control[["aneuploidy_cov"]]

## run the model
aneuploidy_chr18 = runMadSeq(hetero=aneuploidy_hetero, 
                             coverage=aneuploidy_normed_cov, 
                             target_chr="chr18")

## An MadSeq object will be returned
aneuploidy_chr18

The MadSeq object from the runMadSeq function contains:

• The posterior distribution of the selected model
• The delta BIC value between selected model and other models

Note: The value of delta BIC suggests the strength of the confidence of the selected model against other models. In our model, you can set a threshold to get high confidence result, usually it's 20 in our testing cases. We summarize it as follows

BIC = c("[0,10]","(10,20]",">20")
evidence = c("Probably noisy data","Could be positive",
             "High confidence")
table = data.frame(BIC,evidence)
library(knitr)
kable(table,col.names =c("deltaBIC","Evidence against higher BIC") ,align="c")

Visualize the selected model

There are a group of plot functions to plot the output MadSeq object from the runMadSeq.

## plot the posterior distribution for all the parameters in selected model
plotMadSeq(aneuploidy_chr18)

## plot the histogram for the estimated fraction of aneuploidy
plotFraction(aneuploidy_chr18, prob=0.95)

## plot the distribution of AAF as estimated by the model
plotMixture(aneuploidy_chr18)

Description of All Parameters

parameters = c("f","m","mu[1]","mu[2]","mu[3] (LOH model)",
               "mu[3] (meiotic trisomy model)","mu[4]","kappa","p[1]","p[2]",
               "p[3]","p[4]","p[5]","m_cov","p_cov","r_cov")
explains = c("Fraction of mosaic aneuploidy",
             "The midpoint of the alternative allele frequency (AAF) for all heterozygous sites",
             "Mean AAF of mixture 1: the AAFs of this mixture shifted from midpoint to some higher values", 
             "Mean AAF of mixture 2: the AAFs of this mixture shifted from midpoint to some lower values",
             "Mean AAF of mixture 3: In LOH model, mu[3] indicates normal sites without loss of heterozygosity",
             "Mean AAF of mixture 3: In meiotic model, the AAFs of this mixture shifted from 0 to some higher value",
             "Mean AAF of mixture 4: the AAFs of this mixture shifted from 1 to some lower value (only in meiotic model)",
             "Indicate variance of the AAF mixtures: larger kappa means smaller variance",
             "Weight of mixture 1: indicate the proportion of heterozygous sites in the mixture 1",
             "Weight of mixture 2: indicate the proportion of heterozygous sites in the mixture 2",
             "Weight of mixture 3: indicate the proportion of heterozygous sites in the mixture 3 (only in LOH and meiotic model)",
             "Weight of mixture 4: indicate the proportion of heterozygous sites in the mixture 4 (only in meiotic model)",
             "Weight of outlier component: the AAF of 1% sites might not well behaved, so these sites are treated as noise.",
             "Mean coverage of all the sites from the chromosome, estimated from a negative binomial distribution",
             "Prob of the negative binomial distribution for the coverage",
             "Another parameter (r) for the negative binomial disbribution of the coverage, small r means large variance")
table = data.frame(parameters,explains)
kable(table,col.names =c("parameters","description") ,align="c")

MADSEQ documentation built on Nov. 8, 2020, 6:51 p.m.