Build a qseaSet from a sampleTable that lists BAM files and sample
metadata, tile the genome into fixed windows, optionally remove blacklisted
regions, and load coverage and (optionally) CNV information according to the
selected methods. This is typically the first step in a mesa/qsea workflow.
makeQset(
sampleTable,
BSgenome,
chrSelect,
windowSize = 300,
CNVwindowSize = 1e+06,
fragmentType = NULL,
fragmentLength = NULL,
fragmentSD = NULL,
CNVmethod = "HMMdefault",
coverageMethod = "PairedAndR1s",
minMapQual = 10,
minInsertSize = 70,
maxInsertSize = 1000,
minReferenceLength = 30,
badRegions = NULL,
properPairsOnly = FALSE,
hmmCopyGC = NULL,
hmmCopyMap = NULL,
maxPatternDensity = 0.05,
enrichmentMethod = "blind1-15",
parallel = getMesaParallel()
)data.frame.
Must contain sample_name, file_name, and group. If CNVmethod
requires inputs, also input_file. Additional columns are preserved as
sample metadata.
Default: none (must be supplied).
character(1) or NULL.
BSgenome package string (e.g., "BSgenome.Hsapiens.NCBI.GRCh38").
See BSgenome::available.genomes().
Default: NULL (must be supplied).
integer() or character().
Chromosomes to include.
Default: 1:22.
integer(1).
Window size in bp for tiling the genome.
Default: 300.
integer(1).
Copy number variation (CNV) window size (bp).
Default: 1e6.
character(1) or NULL.
If "Sheared" or "cfDNA", sets defaults for fragmentLength/fragmentSD;
otherwise both must be supplied explicitly.
Default: NULL.
numeric(1) or NULL.
Average fragment length (bp). This is used to extend single-end reads if
used, but also strongly affects the CpG density calculation, which will
strongly impact the beta value estimation.
Default: NULL (inferred from fragmentType or must be provided).
numeric(1) or NULL.
Standard deviation of fragment length (bp). Used to calculate the window
based CpG density calculation, which will strongly impact the beta value
estimation.
Default: NULL (inferred from fragmentType or must be provided).
character(1).
One of "HMMdefault", "qseaInput", "MeCap", "None".
Default: "HMMdefault".
character(1).
One of "PairedAndR1s" (mesa) or "qseaPaired" (qsea).
Default: "PairedAndR1s".
integer(1).
Minimum MAPQ to retain a read. For "PairedAndR1s", a pair is kept if
either end passes when properly paired and MAPQ tags are set.
Default: 10.
integer(1).
Minimum absolute insert size for proper pairs (bp).
Only applies to coverageMethod = "PairedAndR1s", and applies to Input
samples as well as MeCap.
Default: 70.
integer(1).
Maximum absolute insert size for proper pairs (bp).
Only applies to coverageMethod = "PairedAndR1s", and applies to Input
samples as well as MeCap.
Default: 1000.
integer(1).
A minimum mapped distance on the genome to keep an read. This is to avoid
very short reads that may be due to mapping artefacts. Note that with
defaultparameters bwa will allow reads where only 19bps have mapped to the
genome.
Only applies to coverageMethod = "PairedAndR1s", and applies to Input
samples as well as MeCap.
Default: 30.
GRanges or NULL.
Genomic regions to exclude (blacklist). These regions will be removed from
the output. Note a set of GRCh38 windows identified by ENCODE are provided
in the ENCODEbadRegions data object.
Default: NULL.
logical(1).
If TRUE, use only proper pairs (stricter size selection).
Default: FALSE.
data.frame or NULL.
GC content per CNV bin (size = CNVwindowSize) for HMMcopy.
Required if using the CNVmethod = "HMMdefault" option unless you are
using hg38/GRCh38 with provided default bin sizes (50kb, 500kb, 1Mb).
Must be a dataframe with columns chr, start, end, and gc.
Default: NULL.
data.frame or NULL.
Mappability per CNV bin (size = CNVwindowSize) for HMMcopy.
Required if using the CNVmethod = "HMMdefault" option unless you are
using hg38/GRCh38 with provided default bin sizes (50kb, 500kb, 1Mb).
Must be a dataframe with columns chr, start, end, and map.
Default: NULL.
numeric(1).
Maximum pattern density (e.g., CpG) to consider a window for background
offset calculation. Regions above this are excluded from the calculation.
This may need to be set higher if you are not considering many windows or
have a high average genomic CG content.
Default: 0.05.
character(1).
Method for calculating enrichment for target pattern (e.g. CpG).
See addNormalisation() for details.
Default: "blind1-15".
logical(1).
Use the registered BiocParallel backend (TRUE) or run serially (FALSE).
See BiocParallel::register().
Default: getMesaParallel().
A qseaSet containing:
tiled (and optionally blacklisted) genomic windows;
the input sampleTable as sample metadata;
loaded coverage according to coverageMethod;
CNV results depending on CNVmethod;
CpG density and enrichment parameters;
recorded processing parameters in @parameters.
Inputs & requirements
sampleTable must include at least:
sample_name (unique per sample),
file_name (path to MeDIP/capture BAM),
group (string label).
If the selected CNVmethod uses inputs ("HMMdefault", "qseaInput"),
also include input_file (path to matched Input BAM).
BAMs should be coordinate-sorted. For paired data, running samtools fixmate
is recommended so MAPQ tags are available (see minMapQual note below).
BSgenome must be an installed genome package string (e.g.
"BSgenome.Hsapiens.NCBI.GRCh38"), which provides chromosome lengths and
sequences for tiling windows.
Windows & blacklist
Chromosomes in chrSelect are tiled into non-overlapping windows of size
windowSize (bp). Windows overlapping badRegions (if provided) are
removed.
Coverage loading
coverageMethod = "PairedAndR1s" loads proper pairs and, optionally,
high-MAPQ R1 singletons; governed by minMapQual, minInsertSize,
maxInsertSize, properPairsOnly, minReferenceLength.
coverageMethod = "qseaPaired" delegates to qsea::addCoverage() with
paired settings.
fragmentType can set defaults for fragmentLength/fragmentSD
("Sheared"
or "cfDNA"); otherwise supply both explicitly.
CNV options
"qseaInput": run qsea::addCNV() using input_file BAMs.
"HMMdefault": use mesa's addHMMcopyCNV() with supplied defaults for
GRCh38
GC/mappability tracks (objects like gc_hg38_1000kb, map_hg38_1000kb,
chosen by CNVwindowSize).
"MeCap": estimate CNV from MeCap BAMs (file_name).
"None": do not compute CNV.
Parallelisation
If parallel = TRUE, computation uses the currently registered
BiocParallel backend (see BiocParallel::register()); otherwise runs
serially.
Use setMesaParallel() to toggle package-wide.
Additional processing
Adds CpG density via qsea::addPatternDensity(pattern = "CG"), and
enrichment parameters, allowing for estimation of beta values.
# A tumour/normal pair of MeDIP-seq BAMs from MEDIPSData, restricted to
# chr22 to keep the example short
if (requireNamespace("MEDIPSData", quietly = TRUE) &&
requireNamespace("BSgenome.Hsapiens.UCSC.hg19", quietly = TRUE)) {
sampleTable <- data.frame(
sample_name = c("Normal1", "Tumour1"),
group = c("Normal", "Tumour"),
file_name = c(
system.file("extdata", "NSCLC_MeDIP_1N_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE),
system.file("extdata", "NSCLC_MeDIP_1T_fst_chr_20_21_22.bam",
package = "MEDIPSData", mustWork = TRUE)
)
)
sampleTable %>%
makeQset(
BSgenome = "BSgenome.Hsapiens.UCSC.hg19",
chrSelect = "chr22",
windowSize = 300,
fragmentLength = 200,
fragmentSD = 50,
CNVmethod = "None",
coverageMethod = "PairedAndR1s",
properPairsOnly = FALSE,
minMapQual = 10
)
}
#> Considering 171015 regions with total size 51304500
#> ==== Creating qsea set ====
#> restricting analysis on chr22
#> Dividing provided regions in 300nt windows
#> Scanning one file at a time. Use e.g. setMesaParallel(nCores = 4) to parallelise this process.
#> Reading NSCLC_MeDIP_1N_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 20294 proper pairs that have neither read with mapq >= 10
#> Dropping 179 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 186, based on 879638 proper pairs.
#> Reading NSCLC_MeDIP_1T_fst_chr_20_21_22.bam, with minimum mapping quality 10.
#> No samtools fixmate MQ (mate quality) tags on the bam file, using R1 MAPQ only.
#> Dropping 12436 proper pairs that have neither read with mapq >= 10
#> Dropping 146 proper pairs that do not lie within [70, 1000]
#> Using high quality R1 unpaired reads as well as pairs
#> Estimating fragment length for 0 unpaired R1 reads as 183, based on 637854 proper pairs.
#> Generated objects
#> Calculated coverage
#> No CNV being calculated
#> Get genomic positions of "CG" ...
#> searching chr22 for "CG"...
#> found 578097 occurances of CG in chr22
#> estimating expected number of CGs per fragment for windows of chr22...
#> ...done
#> selecting windows with low CpG density for background read estimation
#> 2.143% of the windows have enrichment pattern density of at most 0.05 per fragment
#> and are used for background reads estimation
#> qseaSet object generated successfully
#> qseaSet
#> =======================================
#> 2 Samples:
#> Normal1, Tumour1
#> 171015 Regions in 1 chromosomes of BSgenome.Hsapiens.UCSC.hg19