RaMWAS Overview

library("knitr")
# opts_chunk$set(eval = FALSE)
opts_chunk$set(fig.width = 6, fig.height = 6)
library(pander)
panderOptions("digits", 3)

Introduction

RaMWAS provides a complete toolset for methylome-wide association studies (MWAS). It is primarily designed for data from enrichment-based methylation assays, but can be applied to other methylomic data (e.g. Illumina methylation array) as well as other data types like gene expression and genotype data (see vignette).

RaMWAS includes the following major components (steps):

  1. Scanning aligned reads from BAM files.
  2. Calculation of quality control (QC) measures, aggregate QC plots, fragment size distribution estimation.
  3. Creation of CpG score (coverage) matrix, sample normalization.
  4. Principal component analysis.
  5. Methylome-wide association study (MWAS).
  6. Annotation of top findings using the BioMart database.
  7. Multi-marker analysis (methylation risk score) using elastic net.
  8. MWAS with correction for effect of SNPs.

Most steps of RaMWAS are internally parallelized. This is made possible, in part, by the use of the filematrix package for storing the data and accessing it in parallel jobs.

Getting started

Installation

To install the most recent version of RaMWAS please follow instructions at GitHub.com (R 3.3 or newer required). \ To install the Bioconductor version of RaMWAS (R 3.4 or newer requred) use the following commands:

## try http:// if https:// URLs are not supported
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("ramwas")

Loading package and viewing documentation

The package vignettes and reference manual can be viewed online and with the following commands.

library(ramwas) # Loads the package
browseVignettes("ramwas") # Opens vignettes
help(package = "ramwas") # Lists package functions

RaMWAS steps

Pipeline Workflow

To illustrate the main steps of RaMWAS we first create an artificial data set.

suppressPackageStartupMessages(library(ramwas))
# dr = "D:/temp"
library(ramwas)
dr = paste0(tempdir(), "/simulated_project")
ramwas0createArtificialData(
            dir = dr,
            nsamples = 20,
            nreads = 100e3,
            ncpgs = 25e3,
            threads = 2)
cat("Project directory:", dr)

Note. The project directory dr can be set to a more convenient location when running the code.

The function ramwas0createArtificialData created the following files and subdirectories in the project directory.

Each RaMWAS step accepts parameters in the form of a list. Here is the parameter set we will use for all steps below.

param = ramwasParameters(
    dirproject = dr,
    dirbam = "bams",
    filebamlist = "bam_list.txt",
    filecpgset = "Simulated_chromosome.rds",
    cputhreads = 2,
    scoretag = "MAPQ",
    minscore = 4,
    minfragmentsize = 50,
    maxfragmentsize = 250,
    minavgcpgcoverage = 0.3,
    minnonzerosamples = 0.3,
    filecovariates = "covariates.txt",
    modelcovariates = NULL,
    modeloutcome = "age",
    modelPCs = 0,
    toppvthreshold = 1e-5,
    bihost = "grch37.ensembl.org",
    bimart = "ENSEMBL_MART_ENSEMBL",
    bidataset = "hsapiens_gene_ensembl",
    biattributes = c("hgnc_symbol","entrezgene","strand"),
    bifilters = list(with_hgnc_trans_name = TRUE),
    biflank = 0,
    cvnfolds = 5,
    mmalpha = 0,
    mmncpgs = c(5, 10, 50, 100, 500, 1000, 2000)
)

We describe the role of each parameter as they are used below. The complete description of all RaMWAS parameters is available in the parameter vignette.

Scan BAM files and calculate QC indices {#step1}

This step scans all BAM files listed in the filebamlist file, it records read start locations, and calculates a number of quality control metrics. The BAMs must be located in dirbam directory.

Reads are filtered by the scoretag parameter, which is usually either the "MAPQ" field or "AS" tag in the BAM file (see BAM file format). Reads with scores below minscore are excluded.

The minfragmentsize and maxfragmentsize parameters define the minimum and maximum size of DNA fragments that were sequenced. Please note, these parameters are not equal to the read length but instead reflect the length of the DNA fragments that were extracted and sequenced.

The set of CpGs is loaded from the filecpgset file. For more information see the CpG set vignette.

RaMWAS uses cputhreads CPU cores (parallel jobs) to scan BAMs. Hard disk speed can be a bottleneck for this step. If BAMs are stored on a single rotational hard drive, using more than 4 parallel jobs may not provide further speed improvements.

ramwas1scanBams(param)

This creates the following subdirectories in the project directory:

Here is a coverage_by_density plot for the simulated data. It shows higher average CpG score (fragment coverage) for regions with higher CpG densities, up to the saturation point.

pfull = parameterPreprocess(param)
qc = readRDS(paste0(pfull$dirrqc, "/BAM007.qc.rds"))
plot(qc$qc$avg.coverage.by.density)

Note. If a BAM file has previously been scanned, it will not be scanned again in subsequent calls of ramwas1scanBams. This way Step 1 can be rerun multiple times to efficiently include additional BAMs (samples) when they become available.\ Note. The BAM files are no longer needed after this step and the RaMWAS BAM info files are 50 to 100 times smaller than the original BAMs.

Summarize QC measures {#step2}

This step aggregates quality control metrics across all scanned BAM files, produces a number of summary plots and tables, and estimates fragment size distribution.

In practice, multiple BAM files may correspond to the same sample. For simplicity, in the example here each BAM corresponds to one sample with the same name. RaMWAS can be instructed about BAM to sample correspondence via filebam2sample or bam2sample parameters (See parameter vignette).

ramwas2collectqc(param)

The following files and directories are created in the project directory:

After exclusion of BAMs and samples not passing QC, step 2 should be rerun. This ensures that fragment size distribution is estimated using selected data only.

The fragment size distribution estimation plot is presented below. The points indicate the number of reads (y-axis) observed at varying distances from isolated CpGs (x-axis). The red line is the parametic fit for these points. For more information on estimation of fragment size distribution see [@van2013estimation].

qc = readRDS(paste0(pfull$dirqc, "/summary_total/qclist.rds"))
frdata = qc$total$hist.isolated.dist1
estimate = as.double(readLines(
    con = paste0(pfull$dirproject,"/Fragment_size_distribution.txt")))
plotFragmentSizeDistributionEstimate(frdata, estimate)

Calculate CpG score matrix {#step3}

This step creates a CpG score matrix (a.k.a. fragment coverage matrix) for all samples and all CpGs in the CpG set. The samples are defined by either filebam2sample or bam2sample parameter. The CpGs are defined by the filecpgset parameter, see CpG set vignette for more information.

RaMWAS can filter out CpGs with low scores. These CpGs are unmethylated and are unlikely to produce any findings. A CpG is preserved if

For each sample, the CpG scores are affected by the number of sequenced DNA fragments, which varies from sample to sample. To remove this variation, the CpG score matrix is normalized to have the same average score for each sample.

ramwas3normalizedCoverage(param)

This step a new creates directory named coverage_norm_X in the project directory (X is the number of samples, see [Directory names][]) with the following files:

Note. This step created temporary files in the dirtemp directory.

Principal component analysis {#step4}

This step performs principal component analysis (PCA) on the CpG score matrix.

PCA is capturing the main directions of unmeasured sources of variation in the data. The main goal of PCA is estimation of laboratory technical variations which can be used as covariates in the association analyses. PCA can be performed after regressing out known covariates as they represent measured sources of variation, which we need not estimate.

Additionally, large sample loadings of the top PCs can indicate multidimensional outlying samples. Such sample may be excluded from the analysis.

If measured covariates are regressed out prior to conducting the PCA, these covariates are loaded from the file filecovariates. The file can be comma or tab-separated, with sample names in the first column. The artificial data set includes a covariate file covariates.txt.

The parameter modelcovariates names the covariates regressed out before PCA (NULL for none).

ramwas4PCA(param)

This step creates sub-directory PCA_X_cvrts_Y in the directory with score matrix, where X is the number of covariates regressed out and Y is a code distinguishing different sets with the same number of covariates (see [Directory names][]). The sub-directory includes:

The PC plot for artificial data shows one strong component and no outliers in the sample loadings, with first PC clearly capturing sample sex.

eigenvalues = fm.load(paste0(pfull$dirpca, "/eigenvalues"));
eigenvectors = fm.open(
                filenamebase = paste0(pfull$dirpca, "/eigenvectors"),
                readonly = TRUE);
plotPCvalues(eigenvalues)
plotPCvectors(eigenvectors[,1], 1)
plotPCvectors(eigenvectors[,2], 2)
close(eigenvectors)

The file with correlations shows strong correlation of the top PCs with age and sex.

tblcr = read.table(
            file = paste0(pfull$dirpca, "/PC_vs_covs_corr.txt"),
            header = TRUE,
            sep = "\t")
pander(head(tblcr, 3))

The file with p-values indicate statistical significance of these correlations.

Table: Correlations\ in\ PC_vs_covariates_corr.txt\ file.

tblpv = read.table(
            file = paste0(pfull$dirpca, "/PC_vs_covs_pvalue.txt"),
            header = TRUE,
            sep = "\t")
pander(head(tblpv, 3))

Table: P-values\ in\ PC_vs_covariates_pvalue.txt\ file.

Methylome-wide association study (MWAS) {#step5}

This step performs tests for association between normalized CpG scores and the outcome variable named by modeloutcome parameter. The\ analysis corrects for covariates listed in modelcovariates parameter and a number of top principal components (modelPCs). Tests are performed using the linear regression model if the outcome variable is numeric, ANOVA if categorical.

All results are saved in a filematrix. Findings passing the toppvthreshold p-value threshold are recorded in a text file.

ramwas5MWAS(param)

This step creates a sub-directory in the PCA directory named Testing_X_Y_PCs, where X is the name of the outcome variable and Y is the number of top PCs included as covariates (see [Directory names][]). The sub-directory includes:

For the simulated data the QQ-plot for age shows moderate deviation from the diagonal for many CpGs, suggesting many markers with small effects. This is consistent with how the data was generated -- there is weak signal in 1\% of CpGs.

RaMWAS also creates a Manhattan plot with matching Y axis.

mwas = getMWASandLocations(param)
layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(1,2.2))
qqPlotFast(mwas$`p-value`)
man = manPlotPrepare(
        pvalues = mwas$`p-value`,
        chr = mwas$chr,
        pos = mwas$start,
        chrmargins = 0)
manPlotFast(man)

Annotation of top results {#step6}

RaMWAS can annotate top findings using data from biomaRt.

The parameters include:

For detailed description of these parameters and R script for listing allowed values see parameter vignette

Here we annotate top findings with human gene symbols, gene ids, and strand information.

ramwas6annotateTopFindings(param)

The function updates the Top_tests.txt file in the MWAS directory.

chr | position | tstat | pvalue | qvalue | hgnc_symbol | entrezgene | strand :---:|---:|---:|---:|---:|:---:|:---:|:---: chr1 | 15,975,530 | 8.771 | 6.446x10^-8^ | 0.022 | DDI2/DDI2 | 84301/6248 | 1/1 chr1 | 15,975,533 | 8.568 | 9.097x10^-8^ | 0.022 | DDI2/DDI2 | 84301/6248 | 1/1 chr1 | 15,418,248 | -7.654 | 4.571x10^-7^ | 0.071 | KAZN | 23254 | 1

Table: Example\ of\ the\ Top_tests.txt\ file.

Methylation risk score {#step7}

While it is important to find individual CpGs associated with a phenotype, we can often achieve better predictive power by combining information from multiple CpGs.

We take an approach similar to the one proposed by [@horvath2013dna] for predicting biological age from methylation data. In particular, we combine signal across multiple CpGs using the elastic net model [@tibshirani2012strong]. To avoid overfitting and correctly estimate the predictive power of the model we use k-fold cross-validation. Within the cross-validation procedure, for each training set of samples we perform MWAS, select top MWAS sites, train the elastic net, and make predictions for the test samples. The set of predictions is recorded as the MRS.

RaMWAS predicts the outcome variable (modeloutcomes parameter) using top mmncpgs CpGs from the MWAS on the training set of samples. If mmncpgs is a vector of several values, each number of CpGs is tested separately.

The elastic net parameter alpha can be set via mmalpha parameter. The number of folds cvnfolds in the K-fold cross validation is 10 by default.

ramwas7riskScoreCV(param)

This step creates CV_X_folds sub-directory in MWAS directory, where X is the number of folds in k-fold cross validation (see [Directory names][]). It contains:

For the simulated data we get moderately good prediction using just 50 top CpGs.

cv = readRDS(paste0(pfull$dircv, "/rds/CpGs=000050_alpha=0.000000.rds"))
plotPrediction(
        param = pfull,
        outcome = cv$outcome,
        forecast = cv$forecast,
        cpgs2use = 50,
        main = "Prediction success (EN on coverage)")

The correlation of prediction with actual age is maximal when we use 100 top CpGs in elastic net.

cl = readRDS(sprintf("%s/rds/cor_data_alpha=%f.rds",
                    pfull$dircv,
                    pfull$mmalpha))
plotCVcors(cl, pfull)

CpG-SNP analysis

Point mutated CpG sites, called CpG-SNPs [@shoemaker2010allele], are interesting sites for human diseases because in addition to the sequence variation they may show individual differences in DNA methylation. RaMWAS can perform CpG-SNP analyses if SNP data from the same subjects/samples is also available. These tests are performed using a regression model that assesses whether the case-control differences are proportional to the number of CpG sites [@van2015whole].

This type of analysis is explained in the CpG-SNP vignette.

Directory names

The names of CpG score/PCA/MWAS/MRS directories can be recovered by first calling parameterPreprocess function.

pfull = parameterPreprocess(param)
cat("CpG score directory:", "\n", pfull$dircoveragenorm,"\n")
cat("PCA directory:", "\n", pfull$dirpca, "\n")
cat("MWAS directory:", "\n", pfull$dirmwas, "\n")
cat("MRS directory:", "\n", pfull$dircv, "\n")
cat("CpG-SNP directory:", "\n", pfull$dirSNPs, "\n")

Version information

sessionInfo()

References



Try the ramwas package in your browser

Any scripts or data that you put into this service are public.

ramwas documentation built on Oct. 31, 2019, 2:11 a.m.