Spiky: Analysing cfMeDIP-seq data with spike-in controls

Samantha L Wilson and Lauren M Harmon

February 8, 2021

Introduction

To meet the need for a reference control in cell-free methylated DNA immunoprecipitation-sequencing (cfMeDIP-seq)1,2 experiments, we designed spike-in controls and ligated unique molecular indexes (UMI) to adjust for PCR bias, and immunoprecipitation bias caused by the fragment length, G+C content, and CpG density of the DNA fragments that are immunoprecipitated3. This enables absolute quantification of methylated DNA in picomoles, while retaining epigenomic information that allows for sensitive, tissue-specific detection as well as comparable results between different experiments. We designed DNA fragments with 2x3x3x3=54 combinations of methylation status (methylated and unmethylated), fragment length in basepair (bp) (80 bp, 160 bp, 320 bp), G+C content (35%, 50%, 65%), and fraction of CpGs within a fragment (1 CpG/ 80 bp, 1 CpG/ 40 bp, 1 CpG/ 20 bp). Spiky was developed for analyzing DNA methylation of cell-free DNA obtained from cfMeDIP-seq method using reference ‘spike-in’ controls. This package will:

Installation

Install and load the spiky package from Bioconductor.

#To install this package, start R (version "3.6" or later) and enter:
  #if (!requireNamespace("BiocManager", quietly = TRUE))
  #  install.packages("BiocManager")
  #
  #BiocManager::install("spiky")

library(spiky)

Load spike database, or create your own with process_spikes().

Input: A Fasta file, GRanges, or dataframe of spike-in sequences, and a vector of booleans (0 or 1) describing whether each spike-in sequence is methylated.

Output: The output contains a DataFrame with the following columns:

sequence (DNAStringSet) methylated (boolean) CpGs (integer) fmol (numeric) molmass (numeric) GCfrac (numeric) OECpG (numeric) conc (numeric) *frag_grp (character)

If you are using the same synthetic spike-in sequences as described in the manuscript, you may load the spike-in sequence database using:

data(spike)

If you are using custom spike-ins, you can create your own spike-in sequence database using the process_spikes() function, which accepts as input a Fasta file, GRanges, or dataframe, and a vector of booleans (0 or 1) describing whether each spike-in sequence is methylated. The input Fasta file can also be generated from BAM header info as follows:

sb <- system.file("extdata", "example.spike.bam", package="spiky",              mustWork=TRUE)
outFasta <- paste(system.file("extdata", package="spiky", mustWork=TRUE),"/spike_contigs.fa",sep="")
show(generate_spike_fasta(sb, spike=spike,fa=outFasta))
#> Wrote /tmp/RtmpV4lj6p/Rbuild3d8bb64a883e32/spiky/inst/extdata/spike_contigs.fa
#> DNAStringSet object of length 52:
#>      width seq                                              names               
#>  [1]   160 CTTTACTACTGAATGTAAGCTCT...TTTACTATAACCGATTACACAT 160b_2C_35G-2_meth
#>  [2]   160 GTAACATGGTTACCACTGGGACC...ACTAGCCGTGTCCCAACCTCAT 160b_2C_50G-2_meth
#>  [3]   160 GACTCCTCCCTAGGCCCCCATGG...CATGCAGGCCCCCCGCTCCATC 160b_2C_65G-2_meth
#>  [4]   160 GTATAATCATAACAAAGGCCTAA...AGCGACTAAGATCTCAGAATTA 160b_4C_35G-2_mod...
#>  [5]   160 AGCCTTGGACGTGAGTCTCTGTT...ACAATTGTCAGGGCCCTCCAGT 160b_4C_50G-1_meth
#>  ...   ... ...
#> [48]    80 ACAACACCCTCCACCCAATACTT...AAGTCAGTCAAATGCCTGTAAC 80b_2C_50G-1
#> [49]    80 CACCTTGAGACCTCCAGAGGGGG...GTGCGCAGGGGGGAAGGGGGGC 80b_2C_65G-1
#> [50]    80 AAGGCATTACTTATCTAATCAAT...TTTGTACTCGTAGACGAAATTG 80b_4C_35G-1
#> [51]    80 AGTCATCAGCATATTGTCAGTAC...TACTCCTAGTGGGCTGCGTGGT 80b_4C_50G-1
#> [52]    80 TTGGGAGGCTCTGGACTGGGGCA...CGTCCCCCCCCATCCTCTCCGC 80b_4C_65G-1

The spike-in database can then be created with this input Fasta.

spikes <- system.file("extdata", "spikes.fa", package="spiky", mustWork=TRUE)
spikemeth <- spike$methylated
process_spikes(spikes, spikemeth)
#> DataFrame with 52 rows and 5 columns
#>                               sequence methylated      CpGs    GCfrac     OECpG
#>                         <DNAStringSet>  <numeric> <integer> <numeric> <numeric>
#> 80b_1C_35G-1   TAGGATATAG...AAATTATGCT          0         1      0.35      0.29
#> 80b_1C_35G-2   TGTCTAAATT...AAGAACATAT          1         1      0.35      0.29
#> 80b_2C_35G-1   TCTAATACTC...AATCCATAAG          0         2      0.35      0.57
#> 80b_2C_35G-2   CTCAAATATA...CAATAACACT          1         2      0.35      0.57
#> 80b_4C_35G-1   AAGGCATTAC...GACGAAATTG          0         4      0.35      1.14
#> ...                                ...        ...       ...       ...       ...
#> 320b_4C_65G-2  GCAATTGATG...AAGAAGCTAA          0         4      0.35      0.29
#> 320b_8C_65G-1  CCCATGCATC...TCTTACCAGT          1         8      0.50      0.40
#> 320b_8C_65G-2  GAACTTCCAA...ATGCTATGCA          0         8      0.50      0.40
#> 320b_16C_65G-1 TTCGGCACTT...GCTTAAGAAA          1        16      0.50      0.80
#> 320b_16C_65G-2 TTGGGCCGCC...CTGAGATTGT          0        16      0.50      0.80

Process the input files

Spiky supports input files in either BAM or BEDPE format.

BAM Input

BAM file in standard format (For full details about the BAM format, see https://samtools.github.io/hts-specs/SAMv1.pdf). The BAM must also have an accompanying index file, which can be created using samtools index ${filename.bam}. (http://www.htslib.org/doc/samtools-index.html)

BAM required columns

Output: The output objects will be used downstream in the analysis, including

genomic_bam_path <- system.file("extdata", "example_chr21.bam", package="spiky", mustWork=TRUE)
genomic_coverage <- scan_genomic_contigs(genomic_bam_path,spike=spike)
#> Warning in find_spike_contigs(si, spike = spike): Cannot resolve contig name mismatches with spike database names. Your file may not have any spikes, or the spike database might not contain them.
#> Tiling 300bp bins across the genome...Done.
#> Binning genomic coverage...
#> Binning genomic coverage...Done.
#> Done.
spike_bam_path <- system.file("extdata", "example.spike.bam", package="spiky", mustWork=TRUE)
spikes_coverage <- scan_spike_contigs(spike_bam_path,spike=spike)
#> Summarizing spike-in counts...Done.

BEDPE Input

BEDPE file in standard format. For full details about the BEDPE format, see Bedtools documentation (https://bedtools.readthedocs.io/en/latest/content/general-usage.html#bedpe-format). In short, for a pair of ranges 1 and 2, we have fields chrom1, start1, end1, chrom2, start2, end2, and (optionally) name, score, strand1, strand2, plus any other user defined fields that may be included (these are not yet supported by Spiky). For example, two valid BEDPE lines are:

chr1 100 200 chr5 5000 5100 bedpe_example1 30

chr9 900 5000 chr9 3000 3800 bedpe_example2 99 + -

The BEDPE must also have an accompanying index file, which can be created using Bedtools, as in the example shown below, where ${file} represents the name of a BEDPE file.

bedtools sort -i \({file} > sorted_\){file} bgzip sorted_\({file} tabix sorted_\){file}.gz

Output: The output objects will be used downstream in the analysis, including

genomic_bedpe_path <- system.file("extdata", "example_chr21_bedpe.bed.gz", package="spiky", mustWork=TRUE)
genomic_coverage <- scan_genomic_bedpe(genomic_bedpe_path,genome="hg38")
#> Tiling 300bp bins across the genome...Done.
#> Binning genomic coverage...
#> Binning genomic coverage...Done.
#> Done.
spike_bedpe_path <- system.file("extdata", "example_spike_bedpe.bed.gz", package="spiky", mustWork=TRUE)
spikes_coverage <- scan_spike_bedpe(spike_bedpe_path,spike=spike)
#> Warning in .toPairs(res): NAs introduced by coercion
#> Summarizing spike-in counts...Done.

Methylation specificity

For each combination of parameters, we designed two distinct spike-in sequences. One to be methylated and one to be unmethylated. The allows us to assess non-specific binding of the monoclonal antibody on a sample-by-sample basis. To calculate methylation specificity we take the number of methylated reads divided by the total number of reads. It is our recommendation that if methylation specificity is <0.98, then the sample should be flagged or removed from analysis as the cfMeDIP performed inadequately.

This calculation is done by the ‘methylation_specificity’ function.

Input: The output of the ‘scan_spike_contigs’ or ‘scan_spike_bedpe’ functions

Output: methylation specificity mean and median

Example

##Calculate methylation specificity
methyl_spec <- methylation_specificity(spikes_coverage,spike=spike)
print(methyl_spec)
#> $mean
#> [1] 0.8856087
#> 
#> $median
#> [1] 0.9708432

Fit a Gaussian model to predict the molar amount of DNA sequences

For each batch of samples, the coefficients used in the Gaussian generalized linear model will differ. The ‘model_glm_pmol’ will calculate these coefficients and output the model to be used to calculate molar amount (picomoles) on the user’s DNA sequences of interest. We assume that all DNA sequences of interest are methylated after undergoing cfMeDIP-seq. As such, we build the Gaussian generalized linear model on only the methylated spike-in control fragments. A generated Bland-Altman plot will visualize how well the model performs.

Input: The output of the ‘scan_spike_contigs’ or ‘scan_spike_bedpe’ functions

Output:

Example

## Build the Gaussian generalized linear model on the spike-in control data
gaussian_glm <- model_glm_pmol(spikes_coverage,spike=spike)
#> Converting x (a coverage result?) to a data.frame...
summary(gaussian_glm)
#> 
#> Call:
#> glm(formula = conc ~ read_count + fraglen + GC + CpG_3, data = x)
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  3.914e-03  4.919e-04   7.958 8.96e-08 ***
#> read_count  -1.946e-08  9.913e-09  -1.964    0.063 .  
#> fraglen     -1.149e-05  1.295e-06  -8.873 1.50e-08 ***
#> GC           3.702e-06  7.407e-06   0.500    0.622    
#> CpG_3        1.719e-04  2.894e-04   0.594    0.559    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 1.733019e-07)
#> 
#>     Null deviance: 3.3927e-05  on 25  degrees of freedom
#> Residual deviance: 3.6393e-06  on 21  degrees of freedom
#> AIC: -324.54
#> 
#> Number of Fisher Scoring iterations: 2

Calculating molar amount on DNA sequences of interest

For the samples in which the Gaussian generalized linear model was built, we will calculate the molar amount (picomoles) for each DNA sequence of interest.

Input: The output of the ‘scan_genomic_contigs’ or ‘scan_genomic_bedpe’ functions and the Gaussian generalized linear model

Output: sample_pmol_data

Example

# Predict pmol concentration
# To select a genome other than hg38, use BSgenome::available.packages() to find valid BSgenome name
#library("BSgenome.Hsapiens.UCSC.hg38")
sample_data_pmol <- predict_pmol(gaussian_glm, genomic_coverage,bsgenome="BSgenome.Hsapiens.UCSC.hg38",ret="df")
#> Adjusting for bin-level biases...
#> Attempting to load BSgenome.Hsapiens.UCSC.hg38...
#> OK.
#> Done.
#> Starting pmol prediction...
#> Done.
head(sample_data_pmol,n=1)
#>                         chrom range.start range.end
#> chr21:30201601-30201900 chr21    30201601  30201900
#>                                                                                                                                                                                                                                                                                                                             sequence
#> chr21:30201601-30201900 AGGAAGGGGTTCAGTTTCAGTTTTCTGCATATGGCTAGCCAGTTTTCTCAACACCTTTTATTAAATAGGGAATCATTTCCCCATTGCTTGTTTTTGTCAGGTTTGTCAAAGATGAGATGGTTGTAGACGTGCGGTGTTATTTCTGAGGCCTCTGTTCTGTTCCATTGGTCTATATATCTGTTTTGGTATCAGTATCATGCTGTTTTTGTTACTGTAGCCTTGTAGTATAGTTTGAAGTCAGGTAGTGTGATGCCTCCAGCTTTGTTCTTTTTGTTTAGAATTGTCTTGGCTATACAGGCT
#>                         read_count fraglen        GC    CpG_3    pred_conc
#> chr21:30201601-30201900        1.8     300 0.3933333 1.259921 0.0006848288

Adjusting molar amount to binned genomic windows

For our analyses, we binned the genome into 300 bp non-overlapping windows. We then look overlap between fragments in our data with each of the 300 bp genomic windows. We adjust the molar amount (picomoles) by a multiplier. This multiplier is the proportion of overlap between our fragment and the 300 bp window. This is done for every fragment in our sample.

Input: output dataframe produced from predict_pmol

Output: sample_binned_data

Example

sample_binned_data <- bin_pmol(sample_data_pmol)
head(sample_binned_data,n=1)
#>                         chrom range.start range.end
#> chr21:30201601-30201900 chr21    30201601  30201900
#>                                                                                                                                                                                                                                                                                                                             sequence
#> chr21:30201601-30201900 AGGAAGGGGTTCAGTTTCAGTTTTCTGCATATGGCTAGCCAGTTTTCTCAACACCTTTTATTAAATAGGGAATCATTTCCCCATTGCTTGTTTTTGTCAGGTTTGTCAAAGATGAGATGGTTGTAGACGTGCGGTGTTATTTCTGAGGCCTCTGTTCTGTTCCATTGGTCTATATATCTGTTTTGGTATCAGTATCATGCTGTTTTTGTTACTGTAGCCTTGTAGTATAGTTTGAAGTCAGGTAGTGTGATGCCTCCAGCTTTGTTCTTTTTGTTTAGAATTGTCTTGGCTATACAGGCT
#>                         read_count fraglen        GC    CpG_3    pred_conc
#> chr21:30201601-30201900        1.8     300 0.3933333 1.259921 0.0006848288
#>                         adjusted_pred_con
#> chr21:30201601-30201900       0.001232692

Session Info

sessionInfo()
#> R version 4.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [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       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] BSgenome.Hsapiens.UCSC.hg38_1.4.5 spiky_1.10.0                     
#>  [3] Rsamtools_2.20.0                  Biostrings_2.72.0                
#>  [5] XVector_0.44.0                    GenomicRanges_1.56.0             
#>  [7] GenomeInfoDb_1.40.0               IRanges_2.38.0                   
#>  [9] S4Vectors_0.42.0                  BiocGenerics_0.50.0              
#> [11] devtools_2.4.5                    usethis_2.2.3                    
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7                remotes_2.5.0              
#>   [3] rlang_1.1.3                 magrittr_2.0.3             
#>   [5] distributions3_0.2.1        matrixStats_1.3.0          
#>   [7] compiler_4.4.0              mgcv_1.9-1                 
#>   [9] vctrs_0.6.5                 stringr_1.5.1              
#>  [11] profvis_0.3.8               pkgconfig_2.0.3            
#>  [13] crayon_1.5.2                fastmap_1.1.1              
#>  [15] ellipsis_0.3.2              utf8_1.2.4                 
#>  [17] promises_1.3.0              rmarkdown_2.26             
#>  [19] sessioninfo_1.2.2           UCSC.utils_1.0.0           
#>  [21] purrr_1.0.2                 xfun_0.43                  
#>  [23] zlibbioc_1.50.0             cachem_1.0.8               
#>  [25] jsonlite_1.8.8              later_1.3.2                
#>  [27] DelayedArray_0.30.0         BiocParallel_1.38.0        
#>  [29] parallel_4.4.0              R6_2.5.1                   
#>  [31] bslib_0.7.0                 stringi_1.8.3              
#>  [33] rtracklayer_1.64.0          pkgload_1.3.4              
#>  [35] jquerylib_0.1.4             Rcpp_1.0.12                
#>  [37] SummarizedExperiment_1.34.0 knitr_1.46                 
#>  [39] BlandAltmanLeh_0.3.1        tidyselect_1.2.1           
#>  [41] httpuv_1.6.15               Matrix_1.7-0               
#>  [43] splines_4.4.0               rstudioapi_0.16.0          
#>  [45] abind_1.4-5                 yaml_2.3.8                 
#>  [47] codetools_0.2-20            miniUI_0.1.1.1             
#>  [49] curl_5.2.1                  pkgbuild_1.4.4             
#>  [51] lattice_0.22-6              tibble_3.2.1               
#>  [53] Biobase_2.64.0              shiny_1.8.1.1              
#>  [55] withr_3.0.0                 coda_0.19-4.1              
#>  [57] evaluate_0.23               desc_1.4.3                 
#>  [59] survival_3.6-4              urlchecker_1.0.1           
#>  [61] pillar_1.9.0                MatrixGenerics_1.16.0      
#>  [63] generics_0.1.3              rprojroot_2.0.4            
#>  [65] sp_2.1-4                    RCurl_1.98-1.14            
#>  [67] ggplot2_3.5.1               munsell_0.5.1              
#>  [69] scales_1.3.0                xtable_1.8-4               
#>  [71] glue_1.7.0                  tools_4.4.0                
#>  [73] MBA_0.1-0                   BiocIO_1.14.0              
#>  [75] BSgenome_1.72.0             GenomicAlignments_1.40.0   
#>  [77] fs_1.6.4                    mvtnorm_1.2-4              
#>  [79] XML_3.99-0.16.1             grid_4.4.0                 
#>  [81] colorspace_2.1-0            nlme_3.1-164               
#>  [83] GenomeInfoDbData_1.2.12     restfulr_0.0.15            
#>  [85] Formula_1.2-5               cli_3.6.2                  
#>  [87] fansi_1.0.6                 S4Arrays_1.4.0             
#>  [89] dplyr_1.1.4                 bamlss_1.2-4               
#>  [91] gtable_0.3.5                sass_0.4.9                 
#>  [93] digest_0.6.35               SparseArray_1.4.0          
#>  [95] rjson_0.2.21                htmlwidgets_1.6.4          
#>  [97] memoise_2.0.1               htmltools_0.5.8.1          
#>  [99] lifecycle_1.0.4             httr_1.4.7                 
#> [101] mime_0.12

References

1. Shen, S. Y. et al. Sensitive tumour detection and classification using plasma cell-free dna methylomes. Nature 563, 579–583 (2018).

2. Shen, S. Y., Burgener, J. M., Bratman, S. V. & De Carvalho, D. D. Preparation of cfMeDIP-seq libraries for methylome profiling of plasma cell-free dna. Nature protocols 14, 2749–2780 (2019).

3. Wilson, S. L. et al. Sensitive and reproducible cell-free methylome quantification with synthetic spike-in controls. Cell Reports Methods 100294 (2022) doi:https://doi.org/10.1016/j.crmeth.2022.100294.