| Title: | Tools for reading single molecule genomics data |
|---|---|
| Description: | Provides tools for reading of single-molecule genomics data from different file formats in R. The package includes functions for reading collapsed or read-level data from modBAM, bedMethyl and `modkit extract` tabular files, as well as C-to-T conversion based BAM files, as `SummarizedExperiment` objects. |
| Authors: | Michael Stadler [aut, cre] (ORCID: <https://orcid.org/0000-0002-2269-4934>), Charlotte Soneson [aut] (ORCID: <https://orcid.org/0000-0003-3833-2169>), Panagiotis Papasaikas [aut] (ORCID: <https://orcid.org/0000-0002-1640-7636>), Sebastien Smallwood [aut] (ORCID: <https://orcid.org/0009-0006-8466-2459>), Friedrich Miescher Institute for Biomedical Research [cph] |
| Maintainer: | Michael Stadler <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.9 |
| Built: | 2026-10-26 10:48:03 UTC |
| Source: | https://github.com/fmicompbio/SingleMoleculeGenomicsIO |
Add a read-level assay to a SummarizedExperiment object
addReadLevelAssay( se, assayDF, assayName, replaceExisting = FALSE, verbose = FALSE )addReadLevelAssay( se, assayDF, assayName, replaceExisting = FALSE, verbose = FALSE )
se |
A |
assayDF |
A |
assayName |
A character scalar giving the name of the new assay. |
replaceExisting |
A logical scalar indicating whether to replace an existing assay with the same name or not. |
verbose |
Logical scalar. If |
A SummarizedExperiment object with the
new read-level assay added.
Charlotte Soneson
library(SummarizedExperiment) library(S4Vectors) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) # duplicate the 'mod_prob' assay into a new assay named 'new_assay' se <- addReadLevelAssay(se, assayDF = assay(se, "mod_prob"), assayName = "new_assay") assayNames(se) metadata(se)$readLevelDatalibrary(SummarizedExperiment) library(S4Vectors) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) # duplicate the 'mod_prob' assay into a new assay named 'new_assay' se <- addReadLevelAssay(se, assayDF = assay(se, "mod_prob"), assayName = "new_assay") assayNames(se) metadata(se)$readLevelData
Convenience function to extract sequence context around positions of
interest (the rowRanges of
a RangedSummarizedExperiment) and
add them the the SummarizedExperiment's
row data (rowData(se)$sequenceContext). The extracted sequences
will correspond to the regions defined as
resize(rowRanges(x), width = sequenceContextWidth, fix = "center".
Sequence contexts are extracted using extractSeqContext.
addSeqContext(x, sequenceContextWidth, sequenceReference)addSeqContext(x, sequenceContextWidth, sequenceReference)
x |
|
sequenceContextWidth |
A numeric scalar giving the width of the
sequence context to be extracted from the reference
( |
sequenceReference |
A |
A RangedSummarizedExperiment
object with sequence contexts added as a
DNAStringSet object to
rowData(x)$sequenceContext.
Michael Stadler
# load package library(SummarizedExperiment) # file with sequence in fasta format of length 6957060 reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") # define some regions at the end of the reference sequence se <- SummarizedExperiment( assays = matrix(1:3, ncol=1), rowRanges = GRanges( "chr1", IRanges(start = 6957060 - c(4, 2, 0), width = 1, names = c("a", "b", "c")), strand = "-")) # add sequence context (note the padding with N's) rowRanges(se) se <- addSeqContext(se, 7, reffile) rowRanges(se)# load package library(SummarizedExperiment) # file with sequence in fasta format of length 6957060 reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") # define some regions at the end of the reference sequence se <- SummarizedExperiment( assays = matrix(1:3, ncol=1), rowRanges = GRanges( "chr1", IRanges(start = 6957060 - c(4, 2, 0), width = 1, names = c("a", "b", "c")), strand = "-")) # add sequence context (note the padding with N's) rowRanges(se) se <- addSeqContext(se, 7, reffile) rowRanges(se)
Add segemnts for individual reads to the colData of a
SummarizedExperiment object with read-level data.
annotateReadSegments(se, irlList, name)annotateReadSegments(se, irlList, name)
se |
A |
irlList |
A list of |
name |
A character scalar giving the name of the |
A SummarizedExperiment object corresponding to se,
with annotated segments added to colData(se)[[name]]
library(IRanges) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6935400-6936300", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) segs <- list(s1 = IRangesList( "s1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9" = IRanges( start = 6925834, width = 140), "s1-6cf74134-e550-4c02-bd2b-91385422ee25" = IRanges( start = c(6926000, 6926200), width = 140))) se <- annotateReadSegments(se, segs, "nucl") se$nucl lengths(se$nucl$s1)library(IRanges) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6935400-6936300", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) segs <- list(s1 = IRangesList( "s1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9" = IRanges( start = 6925834, width = 140), "s1-6cf74134-e550-4c02-bd2b-91385422ee25" = IRanges( start = c(6926000, 6926200), width = 140))) se <- annotateReadSegments(se, segs, "nucl") se$nucl lengths(se$nucl$s1)
calcReadStats calculates various per-read summary statistics on
modification probabilities or calls from a
SummarizedExperiment object with genomic
positions in rows and samples in columns. addReadStats adds them to
the colData under name. See details for more information on the
statistics that are calculated.
calcReadStats( se, assayName = "mod_prob", stats = NULL, regions = NULL, sequenceContext = NULL, minNobsPpos = 0, minNobsPread = 0, LowConf = 0.7, LagRange = c(12, 64), EntrControl = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE ) addReadStats(se, ..., name = "QC")calcReadStats( se, assayName = "mod_prob", stats = NULL, regions = NULL, sequenceContext = NULL, minNobsPpos = 0, minNobsPread = 0, LowConf = 0.7, LagRange = c(12, 64), EntrControl = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE ) addReadStats(se, ..., name = "QC")
se |
A |
assayName |
A character scalar specifying the assay of |
stats |
Character vector specifying which statistics to calculate.
When set to |
regions |
A |
sequenceContext |
A character vector with sequence context(s)
to include in the calculations. Only positions that match one of the
provided sequence contexts will be included. Sequence contexts can be
provided using IUPAC redundancy codes. The sequence contexts of modified
bases are obtained from |
minNobsPpos |
A numeric scalar value >=1 indicating the minimum coverage on individual positions for them to be included in the calculations. In high coverage data this is an effective filter for removing spurious modbases, typically the result of erroneous basecalling. |
minNobsPread |
A numeric scalar with the minimum number of observed
modifiable bases per read for it to be included in the calculations.
|
LowConf |
A numeric scalar with the minimum call confidence below which calls are considered "low confidence". |
LagRange |
A numeric vector of two values (minimum and maximum) defining the range of lags for the calculation of autocorrelation and partial autocorrelation (see details section). |
EntrControl |
Optional named list with elements
|
BPPARAM |
A |
verbose |
If |
... |
For |
name |
For |
calcReadStats calculates a collection of location/scatter statistics
and information theoretic/signal-processing metrics for the modification
probability, confidence or modification call value vectors across individual
reads (data in assay assayName). Only bases matching the criteria
given byregions, sequenceContext, minNobsPpos and
minNobsPread are included in the calculations. The values of these
filtering parameters are stored in the attribute of the output.
stats selects the summaries to be calculated. Currently available
values are:
: Mean modification probability across the read.
: Fraction of confidently called bases (either modified or unmodified), that are modified. By default, all bases are considered confidently called.
: Mean call confidence across the read.
: Mean call confidence confined to unmodified bases (modifiable bases with modification probability < 0.5).
: Mean call confidence confined to modified bases (modifiable bases with modification probability >= 0.5).
: Fraction of modifiable bases called with low
confidence (call confidence < LowConf).
: Interquartile range of modification probabilities across the read.
: Standard deviation of modification probabilities across the read.
: Sample entropy of the modification probability signal. Sample entropy is a metric assessing the complexity of one-dimensional physiological signals. The higher the sample entropy, the more irregular, unpredictable and therefore complex the signal. See wikipedia:Sample_entropy for more details.
: Autocorrelation of the modification probability values
for lags in the range LagRange. This range typically covers
the signal of nucleosome periodicity.
: Partial autocorrelation of the modification
probability values for lags in the range LagRange. This range
typically covers the signal of nucleosome periodicity.
: raw Read Noise variance estimated as
, where are
lag-1 differences that may skip up to consecutive NAs.
A floor is applied: .
: Read Signal variance
.
: Read Signal-to-Noise Ratio
.
When SNR-related statistics are requested,
a robust noise floor b0 + b1 * mean(x) is fitted per sample (see NoiseVar).
The fitted coefficients are stored in the result metadata (see below). If
SNR-related statistics are not computed, NA/NA placeholders are stored
for uniformity.
For calcReadStats, a SimpleList object with summary statistics
for the samples (columns) in se.
For addReadStats, a SummarizedExperiment
object with summary statistics added to the name column of
colData.
Panagiotis Papapasaikas, Charlotte Soneson, Michael Stadler
# load example data library(SummarizedExperiment) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) readStats <- calcReadStats(se, BPPARAM = BiocParallel::SerialParam()) readStats$s1 se_withReadStats <- addReadStats(se, name = "QC", BPPARAM = BiocParallel::SerialParam()) se_withReadStats$QC$s1 metadata(se_withReadStats$QC$s1)# load example data library(SummarizedExperiment) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) readStats <- calcReadStats(se, BPPARAM = BiocParallel::SerialParam()) readStats$s1 se_withReadStats <- addReadStats(se, name = "QC", BPPARAM = BiocParallel::SerialParam()) se_withReadStats$QC$s1 metadata(se_withReadStats$QC$s1)
All assays with read-level data must have the same number and order of
the reads, which must also agree with the order in se$QC if that
exists. All assays must have the same column names, which must also
agree with the column names of the object, and the sample column
in the colData.
checkSEValidity(se, verbose = FALSE)checkSEValidity(se, verbose = FALSE)
se |
A |
verbose |
A logical scalar. If |
Silently returns NULL. If the object is not valid, an error
will be raised.
Charlotte Soneson
library(GenomicRanges) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, variantPositions = GPos(seqnames = "chr1", pos = c(6940000, 6940500)), BPPARAM = BiocParallel::SerialParam()) checkSEValidity(se)library(GenomicRanges) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, variantPositions = GPos(seqnames = "chr1", pos = c(6940000, 6940500)), BPPARAM = BiocParallel::SerialParam()) checkSEValidity(se)
Report the number of mapped and unmapped alignments in a sorted and indexed bam file.
countAlignmentsInBam(bamfile)countAlignmentsInBam(bamfile)
bamfile |
Character scalar giving the name of the bam file to count alignments in. The file must have an associated index. |
A named list with elements "mapped" and "unmapped",
where the "mapped" element is a named vector with the number
of alignments per chromosome, and the "unmapped" element is
a numeric scalar giving the number of unmapped reads in the bam file.
countAlignmentsInBam(system.file("extdata/6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO"))countAlignmentsInBam(system.file("extdata/6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO"))
For all pairs of bases with modification calls in a read (pair), tabulate the number of pairs at a given distance with a given modification state.
countMismatchStatePairs( bamfile, bamFormat = "QuasR", regions = ".", sequenceContext = "GCH", nAlnsToSample = 0, seqnamesToSampleFrom = character(0), readBaseUnmod = "T", readBaseMod = "C", windowSize = 200, minMapQ = 0, minAlignedLength = 0, seqinfo = NULL, sequenceReference = NULL, BPPARAM = MulticoreParam(4L), verbose = FALSE )countMismatchStatePairs( bamfile, bamFormat = "QuasR", regions = ".", sequenceContext = "GCH", nAlnsToSample = 0, seqnamesToSampleFrom = character(0), readBaseUnmod = "T", readBaseMod = "C", windowSize = 200, minMapQ = 0, minAlignedLength = 0, seqinfo = NULL, sequenceReference = NULL, BPPARAM = MulticoreParam(4L), verbose = FALSE )
bamfile |
Character scalar giving the path to a |
bamFormat |
A character scalar giving the format of |
regions |
A |
sequenceContext |
A character scalar with an odd number of characters,
specifying the genomic context on the plus strand, for which to report
(mis-)matches. The string may contain IUPAC ambiguity codes, e.g.
|
nAlnsToSample |
A numeric scalar. If non-zero, |
seqnamesToSampleFrom |
A character vector with one or several sequence
names (chromosomes) from which to sample alignments from (only used if
|
readBaseUnmod, readBaseMod
|
Character scalars defining the read bases
that are interpreted as unmodified or modified, respectively, when
aligned to the middle base of |
windowSize |
Numeric scalar giving the maximum window size covering pairs of modified bases to consider in pair-counting mode. A window size of 1 corresponds to a single base, a size of 2 to directly adjacent bases, etc. |
minMapQ |
Numeric scalar giving the minimal mapping quality to include alignments in pair-counting mode. |
minAlignedLength |
Numeric scalar giving the minimal alignment length to include alignments in pair-counting mode. |
seqinfo |
|
sequenceReference |
A |
BPPARAM |
A |
verbose |
A logical scalar. If |
A DataFrame with windowSize rows and five columns.
The metadata slot contains header information from bamfile.
Charlotte Soneson, Michael Stadler
bamfile <- system.file("extdata", "BisSeq_quasr_single.bam", package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") tbl <- countMismatchStatePairs(bamfile = bamfile, bamFormat = "QuasR", regions = "chr1:6940000-6955000", sequenceContext = "C", readBaseUnmod = "T", readBaseMod = "C", windowSize = 200, sequenceReference = reffile, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE) tblbamfile <- system.file("extdata", "BisSeq_quasr_single.bam", package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") tbl <- countMismatchStatePairs(bamfile = bamfile, bamFormat = "QuasR", regions = "chr1:6940000-6955000", sequenceContext = "C", readBaseUnmod = "T", readBaseMod = "C", windowSize = 200, sequenceReference = reffile, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE) tbl
For all pairs of bases with modification calls in a read, tabulate the number of pairs at a given distance with a given modification state.
countStatePairs( bamfile, regions = ".", modbase, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), threshUnmod = 0.5, threshMod = 0.5, windowSize = 200, minMapQ = 0, minAlignedLength = 0, BPPARAM = MulticoreParam(4L), verbose = FALSE )countStatePairs( bamfile, regions = ".", modbase, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), threshUnmod = 0.5, threshMod = 0.5, windowSize = 200, minMapQ = 0, minAlignedLength = 0, BPPARAM = MulticoreParam(4L), verbose = FALSE )
bamfile |
Character scalar with the path to a |
regions |
A |
modbase |
Character scalar defining the modified base. |
nAlnsToSample |
A numeric scalar. If non-zero, |
seqnamesToSampleFrom |
A character vector with one or several sequence
names (chromosomes) from which to sample alignments from (only used if
|
threshUnmod, threshMod
|
Numeric scalars defining how to convert
modification probabilities |
windowSize |
Numeric scalar giving the maximum window size covering pairs of modified bases to consider in pair-counting mode. A window size of 1 corresponds to a single base, a size of 2 to directly adjacent bases, etc. |
minMapQ |
Numeric scalar giving the minimal mapping quality to include alignments in pair-counting mode. |
minAlignedLength |
Numeric scalar giving the minimal alignment length to include alignments in pair-counting mode. |
BPPARAM |
A |
verbose |
A logical scalar. If |
A DataFrame with windowSize rows and five columns.
The metadata slot contains header information from bamfile.
Charlotte Soneson, Michael Stadler
modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") res <- countStatePairs(bamfile = modbamfile, regions = "chr1", modbase = "a", windowSize = 300, BPPARAM = BiocParallel::SerialParam()) resmodbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") res <- countStatePairs(bamfile = modbamfile, regions = "chr1", modbase = "a", windowSize = 300, BPPARAM = BiocParallel::SerialParam()) res
Given a time series signal, potentially with missing/unobserved values, this function estimates the level of noise under the assumption of local continuity/smoothness.
estimateNoise(probs, read_pos, k, min_diffs)estimateNoise(probs, read_pos, k, min_diffs)
probs |
Numeric vector of (observed) time series measurements. |
read_pos |
Integer vector of measurement positions in the time series (same length as probs). |
k |
Integer, maximum gap size tolerated when computing |
min_diffs |
Integer, minimum number of lag-1 differences required to produce an estimate (default -1 = auto). |
where are lag-1 differences that may skip
up to k missing values. This follows from error propagation and the assumption
of low varying x in adjacent measurements.
A vector with the following items:
Mean signal
Total variance
Noise variance estimate
Number of adjacent positions used for the estimate
estimateNoise( c(0.1, 0.25, 0.3, 0.45, 0.5, 0.7, 0.7), c(1L, 2L, 3L ,4L, 6L, 8L, 11L), 2L, 1)estimateNoise( c(0.1, 0.25, 0.3, 0.45, 0.5, 0.7, 0.7), c(1L, 2L, 3L ,4L, 6L, 8L, 11L), 2L, 1)
Given time series variance components (total variance, raw noise variance), compute the final noise, signal, and SNR. For noise estimation use one of:
A raw noise variance e.g the result of estimateNoise
A background noise model based on linear predictors
A flooring rule: final noise = max(raw, background)
The background noise model is a general linear noise baseline:
baseline = sum_i beta[i] * feat[i].
estimateSNR(totalVar, noiseRaw, eps, betas, features, noise_mode = "floor")estimateSNR(totalVar, noiseRaw, eps, betas, features, noise_mode = "floor")
totalVar |
Total time series variance. Typically calculated with |
noiseRaw |
noise variance estimate. Typically calculated with |
eps |
Numeric. Lower bound applied to both noise and signal. |
betas, features
|
Coefficients ( |
noise_mode |
Character scalar: one of
|
Named numeric vector with elements: - snr (log2 scale) - signal - noise final estimate used for signal and snr calculation - noise baseline (sum(betas * features); NA if not applied) - noise raw
estimateSNR(totalVar = 1.0, noiseRaw = 1.0, eps = 0.01, betas = c(0.5, 2.0), features = c(2.0, 4.5), noise_mode = "model")estimateSNR(totalVar = 1.0, noiseRaw = 1.0, eps = 0.01, betas = c(0.5, 2.0), features = c(2.0, 4.5), noise_mode = "model")
RangedSummarizedExperiment to single base resolution.Expand the rows of a RangedSummarizedExperiment to single base resolution.
expandSEToBaseSpace( se, regions = NULL, flank = 0L, seqinfo = NULL, keepAssays = getReadLevelAssayNames(se), ignore.strand = TRUE, assayNameAmbig = "mod_prob", verbose = FALSE )expandSEToBaseSpace( se, regions = NULL, flank = 0L, seqinfo = NULL, keepAssays = getReadLevelAssayNames(se), ignore.strand = TRUE, assayNameAmbig = "mod_prob", verbose = FALSE )
se |
|
regions |
A |
flank |
Integer scalar. If |
seqinfo |
|
keepAssays |
Character vector indicating which (read-level) assays to expand to base space. Only these assays will be present in the returned object. |
ignore.strand |
A logical scalar defining whether to ignore the strand
of the |
assayNameAmbig |
A character scalar indicating the assay to use to
decide which row to retain if multiple rows represent the same genomic
position (on different strands). The row with the largest row sum in
this assay is retained. Only used if |
verbose |
Logical scalar indicating whether to print out progress messages. |
A single-base resolution
RangedSummarizedExperiment
corresponding to se. If ignore.strand is TRUE,
the return object is unstranded. If ignore.strand is FALSE,
all positions present in the input object will retain the assigned
strand, while all other positions will be unstranded.
Charlotte Soneson, Michael Stadler
modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se_exp <- expandSEToBaseSpace(se) dim(se) dim(se_exp)modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se_exp <- expandSEToBaseSpace(se) dim(se) dim(se_exp)
Extract one or more (possibly nested) columns from the colData of a SummarizedExperiment object
extractColDataColumns( se, colNames, alignWith = ifelse(any(grepl(":", colNames)), "read", "sample") )extractColDataColumns( se, colNames, alignWith = ifelse(any(grepl(":", colNames)), "read", "sample") )
se |
A |
colNames |
A character vector corresponding to the names of annotation
( |
alignWith |
A character scalar indicating whether the rows in the
returned |
A data.frame with columns sample, read (if
alignWith = "read") and make.names(colNames).
Michael Stadler, Charlotte Soneson
extractfiles <- system.file("extdata", c("modkit_extract_rc_6mA_1.tsv.gz", "modkit_extract_rc_6mA_2.tsv.gz"), package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(extractfiles, modbase = "a", filter = "modkit", BPPARAM = BiocParallel::SerialParam()) df <- extractColDataColumns(se = se, colNames = c("modbase", "readInfo:ref_strand"), alignWith = "read") head(df) dim(df) df <- extractColDataColumns(se = se, colNames = "modbase", alignWith = "sample") head(df) dim(df)extractfiles <- system.file("extdata", c("modkit_extract_rc_6mA_1.tsv.gz", "modkit_extract_rc_6mA_2.tsv.gz"), package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(extractfiles, modbase = "a", filter = "modkit", BPPARAM = BiocParallel::SerialParam()) df <- extractColDataColumns(se = se, colNames = c("modbase", "readInfo:ref_strand"), alignWith = "read") head(df) dim(df) df <- extractColDataColumns(se = se, colNames = "modbase", alignWith = "sample") head(df) dim(df)
This function will extract a sequence context of sequenceContextWidth
bases around the center of the regions defined in x from
sequenceReference.
extractSeqContext(x, sequenceContextWidth, sequenceReference)extractSeqContext(x, sequenceContextWidth, sequenceReference)
x |
A |
sequenceContextWidth |
A numeric scalar giving the width of the
sequence context to be extracted from the reference
( |
sequenceReference |
A |
A DNAStringSet object of the same length
as x with extracted sequence context. All elements are guaranteed
to have identical length (if the sequence context extends to before the
start or beyond the end of a reference sequence, it will be padded with
'N' bases.
Michael Stadler
# file with sequence in fasta format of length 6957060 reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") # define some regions at the end of the reference sequence regions <- GenomicRanges::GRanges( "chr1", IRanges::IRanges(start = 6957060 - c(4, 2, 0), width = 1, names = c("a", "b", "c"))) # extract sequence context (note the padding with N's) extractSeqContext(regions, 7, reffile)# file with sequence in fasta format of length 6957060 reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") # define some regions at the end of the reference sequence regions <- GenomicRanges::GRanges( "chr1", IRanges::IRanges(start = 6957060 - c(4, 2, 0), width = 1, names = c("a", "b", "c"))) # extract sequence context (note the padding with N's) extractSeqContext(regions, 7, reffile)
Filter positions based on any combination of sequence context,
coverage, repetition of the same position, and presence of non-NA values.
Filters are applied in the order they are specified to the filters
argument. Any filter type can be repeated an arbitrary number of times.
filterPositions( se, filters = c("sequenceContext", "coverage", "all.na"), sequenceContext = NULL, assayNameCov = "Nvalid", minCov = 1, minNbrSamples = NULL, regions = NULL, seqinfo = NULL, assayNameAmbig = "Nvalid", assayNameNA = "mod_prob" )filterPositions( se, filters = c("sequenceContext", "coverage", "all.na"), sequenceContext = NULL, assayNameCov = "Nvalid", minCov = 1, minNbrSamples = NULL, regions = NULL, seqinfo = NULL, assayNameAmbig = "Nvalid", assayNameNA = "mod_prob" )
se |
A |
filters |
A character vector. All values must be one of
|
sequenceContext |
A character vector with sequence contexts to
retain. To apply this filter, the |
assayNameCov |
A character scalar indicating the assay to use to
define the coverage. If this is a read-level assay, coverage is first
calculated using |
minCov |
A numeric scalar indicating the lowest acceptable coverage in order to keep a position. |
minNbrSamples |
A numeric scalar, or |
regions |
A |
seqinfo |
|
assayNameAmbig |
A character scalar indicating the assay to use to decide which row to retain if multiple rows represent the same genomic position (on different strands). The row with the largest row sum in this assay is retained. |
assayNameNA |
A character scalar indicating the assay to use as the basis for filtering out positions with NA values across all reads. This should be a read level assay. |
A filtered SummarizedExperiment.
Charlotte Soneson
modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6920000-6940000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- flattenReadLevelAssay(se) se <- addSeqContext(se, sequenceContextWidth = 3, sequenceReference = reffile) sefilt <- filterPositions(se, c("sequenceContext", "coverage", "all.na"), minCov = 5, sequenceContext = "TAG")modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6920000-6940000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- flattenReadLevelAssay(se) se <- addSeqContext(se, sequenceContextWidth = 3, sequenceReference = reffile) sefilt <- filterPositions(se, c("sequenceContext", "coverage", "all.na"), minCov = 5, sequenceContext = "TAG")
Filter reads
filterReads( se, assayName = "mod_prob", readInfoCol = "readInfo", qcCol = "QC", minQscore = 0, maxEntropy = Inf, minSNR = -Inf, maxFracLowConf = 1, minReadLength = 0, minAlignedLength = 0, minAlignedFraction = 0, minCoveredFraction = 0, region = NULL, prune = TRUE, onlyStats = FALSE, removeAllNApos = TRUE )filterReads( se, assayName = "mod_prob", readInfoCol = "readInfo", qcCol = "QC", minQscore = 0, maxEntropy = Inf, minSNR = -Inf, maxFracLowConf = 1, minReadLength = 0, minAlignedLength = 0, minAlignedFraction = 0, minCoveredFraction = 0, region = NULL, prune = TRUE, onlyStats = FALSE, removeAllNApos = TRUE )
se |
A |
assayName |
A character scalar providing the name of a read-level
assay in |
readInfoCol |
A character scalar providing the name of the column in
|
qcCol |
A character scalar providing the name of the column in
|
minQscore |
A numeric scalar representing the smallest acceptable read-level Qscore. Reads with Qscore below this value will be filtered out. |
maxEntropy |
A numeric scalar representing the largest acceptable read-level entropy. Reads without modified-base calls or with entropy above this value will be filtered out. |
minSNR |
A numeric scalar representing the smallest acceptable read Signal to Noise ratio (SNR). Reads with SNR below this value will be filtered out. |
maxFracLowConf |
A numeric scalar representing the maximally acceptable fraction of low-confidence modified base calls in a read. Reads without modified-base calls or with a fraction of low confidence calls greater than this value will be filtered out. |
minReadLength |
A numeric scalar representing the smallest acceptable read length. Reads that are shorter than this value will be filtered out. |
minAlignedLength |
A numeric scalar representing the smallest acceptable aligned length. Reads with aligned length shorter than this value will be filtered out. |
minAlignedFraction |
A numeric scalar representing the smallest acceptable aligned fraction of a read. Reads where the aligned fraction is smaller than this value will be filtered out. |
minCoveredFraction |
A numeric scalar giving the minimal fraction of
|
region |
A |
prune |
A logical scalar. If |
onlyStats |
A logical scalar. If |
removeAllNApos |
A logical scalar. If |
If onlyStats is FALSE, a filtered
SummarizedExperiment object. The metadata of this
object contains a slot named filteredOutReads, which tabulate all
reads that are filtered out, together with the reason(s) for exclusion.
If onlyStats is TRUE, only this table is returned.
Charlotte Soneson, Michael Stadler
library(SummarizedExperiment) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfile, regions = "chr1:6920000-6995000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, name = "QC", BPPARAM = BiocParallel::SerialParam()) ## Filter se ## ... by quality score and aligned length sefilt <- filterReads(se, minQscore = 14, minAlignedLength = 10000) ## ... by covered fraction of the specified region sefilt <- filterReads(se, minCoveredFraction = 1, region = "chr1:6930000-6935000") ## Only calculate filter stats filtstats <- filterReads(se, minQscore = 14, minAlignedLength = 10000, onlyStats = TRUE) filtstats ## Visualize filter stats in UpSet plot, e.g. with ComplexUpset if (require(SimpleUpset)) { SimpleUpset::simpleUpSet(as.data.frame(filtstats$s1)) }library(SummarizedExperiment) modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfile, regions = "chr1:6920000-6995000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, name = "QC", BPPARAM = BiocParallel::SerialParam()) ## Filter se ## ... by quality score and aligned length sefilt <- filterReads(se, minQscore = 14, minAlignedLength = 10000) ## ... by covered fraction of the specified region sefilt <- filterReads(se, minCoveredFraction = 1, region = "chr1:6930000-6935000") ## Only calculate filter stats filtstats <- filterReads(se, minQscore = 14, minAlignedLength = 10000, onlyStats = TRUE) filtstats ## Visualize filter stats in UpSet plot, e.g. with ComplexUpset if (require(SimpleUpset)) { SimpleUpset::simpleUpSet(as.data.frame(filtstats$s1)) }
infile to outfile if they pass filter
criteria.For each bam file in infiles, parse alignments, calculate read
statistics and write the alignment to the corresponding output file from
outfiles if the read passes all criteria defined by the filtering
arguments. Filters are processed hierarchically: If a read does not pass a
given filter, the remaining filters will not be examined and the processing
continues with the next read.
The filters are examined in this order: keepUnmapped,
keepSecondary, keepSupplementary, minReadLength,
minAlignedLength, minAlignedFraction, minQscore,
minSNR, maxFracLowConf, maxEntropy.
filterReadsModBam( infiles, outfiles, modbase, indexOutfiles = TRUE, overwriteOutfiles = FALSE, keepUnmapped = TRUE, keepSecondary = TRUE, keepSupplementary = TRUE, minReadLength = 0, minAlignedLength = 0, minAlignedFraction = 0, minQscore = 0, minSNR = -Inf, maxFracLowConf = 1, maxEntropy = Inf, noiseCoef = c(NA_real_, NA_real_), LowConf = 0.7, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )filterReadsModBam( infiles, outfiles, modbase, indexOutfiles = TRUE, overwriteOutfiles = FALSE, keepUnmapped = TRUE, keepSecondary = TRUE, keepSupplementary = TRUE, minReadLength = 0, minAlignedLength = 0, minAlignedFraction = 0, minQscore = 0, minSNR = -Inf, maxFracLowConf = 1, maxEntropy = Inf, noiseCoef = c(NA_real_, NA_real_), LowConf = 0.7, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )
infiles |
Character vector with name(s) of the input bam file(s). |
outfiles |
Character vector with name(s) of the output bam file(s).
Needs to have the same length as |
modbase |
Character scalar defining the modified base to analyze
(used by |
indexOutfiles |
Logical scalar. If |
overwriteOutfiles |
Logical scalar. If |
keepUnmapped, keepSecondary, keepSupplementary
|
Logical scalars indicating whether to keep unmapped, secondary or supplementary alignments. |
minReadLength |
A numeric scalar representing the smallest acceptable read length. Reads that are shorter than this value will be filtered out. |
minAlignedLength |
A numeric scalar representing the smallest acceptable aligned length. Reads with aligned length shorter than this value will be filtered out. |
minAlignedFraction |
A numeric scalar representing the smallest acceptable aligned fraction of a read. Reads where the aligned fraction is smaller than this value will be filtered out. |
minQscore |
A numeric scalar representing the smallest acceptable read-level Qscore. Reads with Qscore below this value will be filtered out. |
minSNR |
Numeric scalar. Minimum acceptable read Signal-to-Noise Ratio (SNR).
Reads with SNR below this value are filtered out. Set to |
maxFracLowConf |
A numeric scalar representing the maximally acceptable fraction of low-confidence modified base calls in a read. Reads without modified-base calls or with a fraction of low confidence calls greater than this value will be filtered out. |
maxEntropy |
A numeric scalar representing the largest acceptable
read-level entropy. Reads without modified-base calls or with entropy
above this value will be filtered out. A value of |
noiseCoef |
Numeric vector of length 2 giving the background noise model
coefficients |
LowConf |
A numeric scalar with the minimum call confidence below which calls are considered "low confidence". |
BPPARAM |
A |
verbose |
Logical scalar. If |
filterReadsModBam is called for its side effect of generating
new modBam files containing the subset of bam records from input modBam
files that pass all filtering criteria. In addition, it returns a
data.frame with one row per infiles giving the numbers of
bam records that were read in total, that were retained in
the outfiles and that were filtered-out by reason of exclusion.
Michael Stadler
modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") filtbamfiles <- tempfile(fileext = rep(".bam", length(modbamfiles))) res <- filterReadsModBam(infiles = modbamfiles, outfiles = filtbamfiles, modbase = "a", indexOutfiles = FALSE, minReadLength = 6746, minAlignedLength = 6896, minAlignedFraction = 0.56, minSNR=-0.768, minQscore = 9.7, maxFracLowConf = 0.11, maxEntropy = 0.29, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE) res unlink(filtbamfiles)modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") filtbamfiles <- tempfile(fileext = rep(".bam", length(modbamfiles))) res <- filterReadsModBam(infiles = modbamfiles, outfiles = filtbamfiles, modbase = "a", indexOutfiles = FALSE, minReadLength = 6746, minAlignedLength = 6896, minAlignedFraction = 0.56, minSNR=-0.768, minQscore = 9.7, maxFracLowConf = 0.11, maxEntropy = 0.29, BPPARAM = BiocParallel::SerialParam(), verbose = TRUE) res unlink(filtbamfiles)
This function will take a SummarizedExperiment
object with read-level footprinting data (e.g. returned by
readModkitExtract or readModBam) and summarize
reads in each sample, for instance to generate modified and total counts at
each position for each sample.
flattenReadLevelAssay( se, assayName = "mod_prob", statistics = c("Nmod", "Nvalid", "FracMod"), modProbThreshold = 0.5, keepReads = TRUE, replaceExisting = TRUE, verbose = FALSE )flattenReadLevelAssay( se, assayName = "mod_prob", statistics = c("Nmod", "Nvalid", "FracMod"), modProbThreshold = 0.5, keepReads = TRUE, replaceExisting = TRUE, verbose = FALSE )
se |
|
assayName |
A character scalar specifying the assay of |
statistics |
Character vector specifying the type of statistics to be
computed. Currently supported values are "Nmod" (number of values per
row in the |
modProbThreshold |
A numeric scalar, indicating the modification probability threshold to use to classify a base as 'modified' or 'unmodified'. |
keepReads |
A logical scalar. If |
replaceExisting |
A logical scalar. If |
verbose |
If |
A SummarizedExperiment object
with the same dimensions as se (positions in rows and samples in
columns), and added assays corresponding to the requested statistics.
Charlotte Soneson, Michael Stadler
SummarizedExperiment for the
returned object type, readModkitExtract for the function
used to read the input files.
exfile <- system.file("extdata", "modkit_extract_rc_6mA_1.tsv.gz", package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(exfile, modbase = "a", BPPARAM = BiocParallel::SerialParam()) se se_summary <- flattenReadLevelAssay(se) se_summaryexfile <- system.file("extdata", "modkit_extract_rc_6mA_1.tsv.gz", package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(exfile, modbase = "a", BPPARAM = BiocParallel::SerialParam()) se se_summary <- flattenReadLevelAssay(se) se_summary
Get base coverage histogram for a BAM file
getBaseCoverageForBam( bamfile, regions = NULL, maxDepth = 200L, method = "full", nThreads = 3L )getBaseCoverageForBam( bamfile, regions = NULL, maxDepth = 200L, method = "full", nThreads = 3L )
bamfile |
A character scalar with the bam file name (and path). |
regions |
Character vector specifying the region(s) for which
to calculate coverage. Each region uses the grammar supported by the
|
maxDepth |
An integer scalar defining the maximal depth to consider. |
method |
Character scalar with the method used for coverage
calculation. |
nThreads |
A numeric scalar with the number of threads used for decompressing BAM records. |
Secondary and supplementary alignments and unmapped reads are not included.
A named numeric vector of length maxDepth + 1, with values at
index i giving the number of positions that were overlapped by
exactly i-1 alignments. Positions overlapped by more than
maxDepth alignments are also added to the value for maxDepth
at index maxDepth + 1.
The algorithm was described in Pedersen BS and Quinlan AR. "Mosdepth: quick coverage calculation for genomes and exomes". Bioinformatics. 2018; 34(5):867-868. https://doi.org/10.1093/bioinformatics/btx699
modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") getBaseCoverageForBam(modbamfile, "chr1", 12L, "full") getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple")modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") getBaseCoverageForBam(modbamfile, "chr1", 12L, "full") getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple")
The names of assays designated as containing read-level data are
extracted from metadata(se)$readLevelData$assayNames.
getReadLevelAssayNames(se)getReadLevelAssayNames(se)
se |
A |
A (possibly empty) character vector with the names of the assays of se containing read-level data.
Charlotte Soneson
modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) getReadLevelAssayNames(se)modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) getReadLevelAssayNames(se)
The names of colData column designated as containing read-level
annotations are extracted from metadata(se)$readLevelData$colDataColumns.
getReadLevelColDataNames(se)getReadLevelColDataNames(se)
se |
A |
A (possibly empty) character vector with the names of the columns of colData(se) containing read-level data.
Charlotte Soneson
modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) getReadLevelColDataNames(se)modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = FALSE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) getReadLevelColDataNames(se)
Get list of read names by sample from SummarizedExperiment object
getReadNamesBySample(se)getReadNamesBySample(se)
se |
A |
A named list with one entry per sample, containing the read names for the respective sample.
Charlotte Soneson
extractfiles <- system.file("extdata", c("modkit_extract_rc_6mA_1.tsv.gz", "modkit_extract_rc_6mA_2.tsv.gz"), package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(extractfiles, modbase = "a", filter = "modkit", BPPARAM = BiocParallel::SerialParam()) getReadNamesBySample(se)extractfiles <- system.file("extdata", c("modkit_extract_rc_6mA_1.tsv.gz", "modkit_extract_rc_6mA_2.tsv.gz"), package = "SingleMoleculeGenomicsIO") se <- readModkitExtract(extractfiles, modbase = "a", filter = "modkit", BPPARAM = BiocParallel::SerialParam()) getReadNamesBySample(se)
This function plots histograms of the various per-read summary statistics
calculated by addReadStats.
plotReadStats( se, readInfoCol = "readInfo", qcCol = "QC", minQscore = NULL, maxEntropy = NULL, maxFracLowConf = NULL, minReadLength = NULL, minAlignedLength = NULL, minAlignedFraction = NULL, minSNR = NULL )plotReadStats( se, readInfoCol = "readInfo", qcCol = "QC", minQscore = NULL, maxEntropy = NULL, maxFracLowConf = NULL, minReadLength = NULL, minAlignedLength = NULL, minAlignedFraction = NULL, minSNR = NULL )
se |
A |
readInfoCol |
A character scalar providing the name of the column in
|
qcCol |
A character scalar providing the name of the column in
|
minQscore, maxEntropy, maxFracLowConf, minReadLength, minAlignedLength, minAlignedFraction, minSNR
|
Numeric scalars representing possible threshold values used in
|
A ggplot object.
Panagiotis Papapasaikas, Charlotte Soneson, Michael Stadler
library(SummarizedExperiment) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) plotReadStats(se)library(SummarizedExperiment) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfile = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam()) se <- addReadStats(se, BPPARAM = BiocParallel::SerialParam()) plotReadStats(se)
bedMethyl fileThis function will read collapsed single-molecule footprinting data
(reads combined per genomic position) from a bedMethyl file.
readBedMethyl( fnames, modbase, nrows = Inf, sampleAnnot = NULL, seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )readBedMethyl( fnames, modbase, nrows = Inf, sampleAnnot = NULL, seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )
fnames |
Character vector with one or several paths of |
modbase |
Character vector defining the modified base for each sample.
If |
nrows |
Only read |
sampleAnnot |
A |
seqinfo |
|
sequenceContextWidth, sequenceReference
|
Define the sequence
context to be extracted around modified bases. By default (
|
BPPARAM |
A |
verbose |
If |
A SummarizedExperiment object
with genomic positions in rows and samples (the unique names of
fnames) in the columns. If sequenceContextWidth != 0,
rowData(x)$sequenceContext will be a DNAStringSet
object with the extracted sequences.
Michael Stadler, Charlotte Soneson
modkit software,
bedMethyl format description,
SummarizedExperiment for the returned object type,
fread for the function used to read the input files,
addSeqContext used to add the sequence context.
bmfile <- system.file("extdata", "modkit_pileup_5mC_1.bed.gz", package = "SingleMoleculeGenomicsIO") readBedMethyl(bmfile, modbase = "m", BPPARAM = BiocParallel::SerialParam())bmfile <- system.file("extdata", "modkit_pileup_5mC_1.bed.gz", package = "SingleMoleculeGenomicsIO") readBedMethyl(bmfile, modbase = "m", BPPARAM = BiocParallel::SerialParam())
Parse read-to-genome mismatches and return a
SummarizedExperiment object with
information on base states. Base mismatches are typically resulting from
single molecule genomics experiments, representing the accessibility of
individual bases: For bisulfite-sequencing experiments, C-to-T mismatches
correspond to unmethylated (inaccessible) bases (readBaseUnmod="T",
readBaseMod="C", see below), while in deaminase-treatment based experiments,
C-to-T mismatches represent modified (accessible) bases
(readBaseUnmod="C", readBaseMod="T"). Secondary and supplementary
alignments are ignored.
readMismatchBam( bamfiles, bamFormat = "QuasR", regions = NULL, sequenceContext = "GCH", readBaseUnmod = "T", readBaseMod = "C", level = "read", overlapAggregation = "maxQscore", sampleAnnot = NULL, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), seqinfo = NULL, sequenceReference = NULL, variantPositions = NULL, trim = FALSE, maxCoverage = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )readMismatchBam( bamfiles, bamFormat = "QuasR", regions = NULL, sequenceContext = "GCH", readBaseUnmod = "T", readBaseMod = "C", level = "read", overlapAggregation = "maxQscore", sampleAnnot = NULL, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), seqinfo = NULL, sequenceReference = NULL, variantPositions = NULL, trim = FALSE, maxCoverage = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )
bamfiles |
Character vector with one or several paths of |
bamFormat |
A character scalar giving the format of |
regions |
A |
sequenceContext |
A character scalar with an odd number of characters,
specifying the genomic context on the plus strand, for which to report
(mis-)matches. The string may contain IUPAC ambiguity codes, e.g.
|
readBaseUnmod, readBaseMod
|
Character scalars defining the read bases
that are interpreted as unmodified or modified, respectively, when
aligned to the middle base of |
level |
Character scalar specifying the level of returned modification data. Supported values are:
|
overlapAggregation |
Character scalar indicating how to aggregate
modification probabilities for positions that are covered by both
mates in a paired-end bam file. Only used if |
sampleAnnot |
A |
nAlnsToSample |
A numeric scalar. If non-zero, |
seqnamesToSampleFrom |
A character vector with one or several sequence
names (chromosomes) from which to sample alignments from (only used if
|
seqinfo |
|
sequenceReference |
A |
variantPositions |
An optional |
trim |
A logical scalar. If |
maxCoverage |
Integer scalar used to increase the maximal coverage for
which samtools pileup will allocate cache memory. Ignored if
|
BPPARAM |
A |
verbose |
Logical scalar. If |
A SummarizedExperiment object
with genomic positions in rows and samples in columns. The assays
depend on the value of the level argument (see above).
Michael Stadler, Charlotte Soneson
bamfile <- system.file("extdata", "BisSeq_quasr_single.bam", package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") se <- readMismatchBam(bamfiles = bamfile, regions = "chr1:6940000-6955000", sequenceReference = reffile, sequenceContext = "NCG", readBaseUnmod = "T", readBaseMod = "C", verbose = TRUE, BPPARAM = BiocParallel::SerialParam())bamfile <- system.file("extdata", "BisSeq_quasr_single.bam", package = "SingleMoleculeGenomicsIO") reffile <- system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO") se <- readMismatchBam(bamfiles = bamfile, regions = "chr1:6940000-6955000", sequenceReference = reffile, sequenceContext = "NCG", readBaseUnmod = "T", readBaseMod = "C", verbose = TRUE, BPPARAM = BiocParallel::SerialParam())
Parse ML and MM tags (see https://samtools.github.io/hts-specs/SAMtags.pdf,
section 1.7) and return a SummarizedExperiment
object with information on modified bases. Implicitly called bases will get
a modification probability of zero. Secondary and supplementary alignments
are ignored.
readModBam( bamfiles, regions = NULL, modbase, level = ifelse(nAlnsToSample == 0 & is.null(variantPositions), "quickread", "read"), sampleAnnot = NULL, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, variantPositions = NULL, modProbThreshold = 0.5, trim = FALSE, maxCoverage = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )readModBam( bamfiles, regions = NULL, modbase, level = ifelse(nAlnsToSample == 0 & is.null(variantPositions), "quickread", "read"), sampleAnnot = NULL, nAlnsToSample = 0, seqnamesToSampleFrom = character(0), seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, variantPositions = NULL, modProbThreshold = 0.5, trim = FALSE, maxCoverage = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )
bamfiles |
Character vector with one or several paths of |
regions |
A |
modbase |
Character vector defining the modified base for each sample.
If |
level |
Character scalar specifying the level of returned modification data. Supported values are:
|
sampleAnnot |
A |
nAlnsToSample |
A numeric scalar. If non-zero, |
seqnamesToSampleFrom |
A character vector with one or several sequence
names (chromosomes) from which to sample alignments from (only used if
|
seqinfo |
|
sequenceContextWidth, sequenceReference
|
Define the sequence
context to be extracted around modified bases. By default (
|
variantPositions |
An optional |
modProbThreshold |
A numeric scalar, indicating the modification
probability threshold to use to classify a base as 'modified' or
'unmodified'. Only used if |
trim |
A logical scalar. If |
maxCoverage |
Integer scalar used to increase the maximal coverage for
which samtools pileup will allocate cache memory. Ignored if
|
BPPARAM |
A |
verbose |
Logical scalar. If |
A SummarizedExperiment object
with genomic positions in rows and samples in columns. The assays
depend on the value of the level argument (see above).
Michael Stadler, Charlotte Soneson
https://samtools.github.io/hts-specs/SAMtags.pdf describing the SAM ML and MM tags for base modifications.
modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam())modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, BPPARAM = BiocParallel::SerialParam())
This function reads one or more modkit extract files and generates a
single SummarizedExperiment object with rows representing positions and
columns representing samples. The SummarizedExperiment object
contains one assay (mod_prob) with modification probabilities for
each position in each read. Unobserved read/position combinations are
represented by NA.
readModkitExtract( fnames, modbase, filter = NULL, nrows = Inf, sampleAnnot = NULL, seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )readModkitExtract( fnames, modbase, filter = NULL, nrows = Inf, sampleAnnot = NULL, seqinfo = NULL, sequenceContextWidth = 0, sequenceReference = NULL, BPPARAM = MulticoreParam(4L, RNGseed = 42L), verbose = FALSE )
fnames |
Character vector with one or several paths of |
modbase |
Character vector defining the modified base for each sample.
If |
filter |
Either |
nrows |
Only read |
sampleAnnot |
A |
seqinfo |
|
sequenceContextWidth, sequenceReference
|
Define the sequence
context to be extracted around modified bases. By default (
|
BPPARAM |
A |
verbose |
If |
A SummarizedExperiment object
with genomic positions in rows and samples in columns. The assay
"mod_prob" contains per-read modification probabilities,
with each column (sample) corresponding to a position-by-read
NaMatrix.
Charlotte Soneson, Michael Stadler
modkit software,
SummarizedExperiment for the returned object type,
fread for the function used to read the input files
extrfile <- system.file("extdata", "modkit_extract_rc_5mC_1.tsv.gz", package = "SingleMoleculeGenomicsIO") ## ... no filtering readModkitExtract(extrfile, modbase = "m", filter = NULL, BPPARAM = BiocParallel::SerialParam()) ## ... modkit filtering readModkitExtract(extrfile, modbase = "m", filter = "modkit", BPPARAM = BiocParallel::SerialParam())extrfile <- system.file("extdata", "modkit_extract_rc_5mC_1.tsv.gz", package = "SingleMoleculeGenomicsIO") ## ... no filtering readModkitExtract(extrfile, modbase = "m", filter = NULL, BPPARAM = BiocParallel::SerialParam()) ## ... modkit filtering readModkitExtract(extrfile, modbase = "m", filter = "modkit", BPPARAM = BiocParallel::SerialParam())
This function takes a supported argument value specifying the
reference sequence (e.g. genome) and returns it as a
DNAStringSet object.
Currently supported argument values are:
: A BSgenome object
: A DNAStringSet object
: A path to a fasta file
refargToDNAStringSet(sequenceReference)refargToDNAStringSet(sequenceReference)
sequenceReference |
A supported value defining the reference sequence. |
A DNAStringSet object.
Michael Stadler, Charlotte Soneson
ref <- refargToDNAStringSet(system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO")) refref <- refargToDNAStringSet(system.file("extdata", "reference.fa.gz", package = "SingleMoleculeGenomicsIO")) ref
GRanges objectThis function takes a character vector with one or several strings
corresponding to genomic regions and converts them to a
GRanges object. The supported forms
of strings are the same as the ones supported by the htslib C API
that is for example used in samtools:
: All of seqname REF
: Seqname REF from START to end of REF
: Seqname REF from 1 to END
: Seqname REF from START to END
: All seqnames from 1 to their ends
Please note that the returned GRanges is generally parallel
to the elements of regions, with the exception of ".":
This special region is only allowed as a length-one regions
argument and will possibly result several returned regions
corresponding to the sequences in seqinfo.
regionStringToGRanges(regions, seqinfo = NULL, maxend = .Machine$integer.max)regionStringToGRanges(regions, seqinfo = NULL, maxend = .Machine$integer.max)
regions |
Character vector with region strings. See details for supported forms. |
seqinfo |
One of |
maxend |
Integer scalar giving the maximal value for END in cases
where it is not given in |
A GRanges object.
Michael Stadler
regionStringToGRanges("chr1:6940000-6955000") regionStringToGRanges("chr1", seqinfo = c(chr1 = 999)) regionStringToGRanges("chr1")regionStringToGRanges("chr1:6940000-6955000") regionStringToGRanges("chr1", seqinfo = c(chr1 = 999)) regionStringToGRanges("chr1")
Regroup reads into new "samples", defined by readGroups.
regroupReads(se, readGroups) regroupReadsByColData(se, colNames, withinSample = FALSE)regroupReads(se, readGroups) regroupReadsByColData(se, colNames, withinSample = FALSE)
se |
A |
readGroups |
A named list, where each element corresponds to a new group of reads (a "sample"). |
colNames |
A character vector corresponding to the names of annotation
( |
withinSample |
A logical scalar, indicating whether the regrouping
should be done within each current sample (column of |
A SummarizedExperiment object
with columns corresponding to the names of readGroups. Only the
read-level assays will be returned (as summary-level assays are likely
no longer consistent with the regrouped read-level assays).
Charlotte Soneson
library(SummarizedExperiment) library(GenomicRanges) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, variantPositions = GPos(seqnames = "chr1", pos = c(6940000, 6940500)), BPPARAM = BiocParallel::SerialParam()) # number of reads per sample lapply(assay(se, "mod_prob"), ncol) # define groups groups <- list(g1 = c("s1-233e48a7-f379-4dcf-9270-958231125563", "s2-d03efe3b-a45b-430b-9cb6-7e5882e4faf8"), g2 = "s1-92e906ae-cddb-4347-a114-bf9137761a8d", g3 = c("s2-034b625e-6230-4f8d-a713-3a32cd96c298", "s1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9")) # regroup reads sere <- regroupReads(se, readGroups = groups) lapply(assay(sere, "mod_prob"), ncol) # regroup by read annotation (variant label) # ... across samples sere <- regroupReadsByColData(se, colNames = "readInfo:variant_label", withinSample = FALSE) sere lapply(assay(sere, "mod_prob"), ncol) # ... within sample sere <- regroupReadsByColData(se, colNames = "readInfo:variant_label", withinSample = TRUE) sere lapply(assay(sere, "mod_prob"), ncol)library(SummarizedExperiment) library(GenomicRanges) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(bamfiles = modbamfiles, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE, variantPositions = GPos(seqnames = "chr1", pos = c(6940000, 6940500)), BPPARAM = BiocParallel::SerialParam()) # number of reads per sample lapply(assay(se, "mod_prob"), ncol) # define groups groups <- list(g1 = c("s1-233e48a7-f379-4dcf-9270-958231125563", "s2-d03efe3b-a45b-430b-9cb6-7e5882e4faf8"), g2 = "s1-92e906ae-cddb-4347-a114-bf9137761a8d", g3 = c("s2-034b625e-6230-4f8d-a713-3a32cd96c298", "s1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9")) # regroup reads sere <- regroupReads(se, readGroups = groups) lapply(assay(sere, "mod_prob"), ncol) # regroup by read annotation (variant label) # ... across samples sere <- regroupReadsByColData(se, colNames = "readInfo:variant_label", withinSample = FALSE) sere lapply(assay(sere, "mod_prob"), ncol) # ... within sample sere <- regroupReadsByColData(se, colNames = "readInfo:variant_label", withinSample = TRUE) sere lapply(assay(sere, "mod_prob"), ncol)
sampleEntropy returns the sample entropy of a time-series signal
see also https://en.wikipedia.org/wiki/Sample_entropy
sampleEntropy(data, m, r, maxStarts = 1000L, nThreads = 1L)sampleEntropy(data, m, r, maxStarts = 1000L, nThreads = 1L)
data |
Numeric vector |
m |
Integer. The pattern (embedding) length: Larger m captures finer structure but sharply reduces the number of matches, so it requires longer data and increases variance. Common choices are m = 2 or 3 for physiological time series. |
r |
Scaling parameter for the filtering factor. The filtering factor is r x standard deviation of the signal |
maxStarts |
Integer specifying the maximum number of signal start
positions to evaluate (default |
nThreads |
Integer giving the number of parallel OpenMP threads to use for calculation. |
This function calculates the Sample Entropy of a time-series vector given as argument. Sample Entropy is used to assess the complexity of physiological time-series signals. This C++ implementation is a modified version of: https://gist.github.com/schochastics/e3684645763e93cbc2ed7d1b70ee5fe6
The Sample Entropy value of the time-series signal
wikipedia:Sample_entropy Multiscale entropy of biological signals
ts <- runif(100, 0, 1) sampleEntropy(ts, m = 2L, r = 0.2)ts <- runif(100, 0, 1) sampleEntropy(ts, m = 2L, r = 0.2)
This function takes read names or indices and subsets them from the
specified assays of RangedSummarizedExperiment
with read-level data. While subsetting of samples (the columns of the
SummarizedExperiment object) can be done easily (e.g. se[, 1]),
the reads are grouped by sample in read-level assays, and this function
provides a more convenient way to subset these nested reads.
subsetReads( se, reads = NULL, randomSubset = NULL, prune = TRUE, invert = FALSE, removeAllNApos = FALSE, assayNameNA = "mod_prob" )subsetReads( se, reads = NULL, randomSubset = NULL, prune = TRUE, invert = FALSE, removeAllNApos = FALSE, assayNameNA = "mod_prob" )
se |
A |
reads |
Defines reads to retain (or remove, if |
randomSubset |
A numeric scalar, enabling sampling of a random set of
reads for each sample. If |
prune |
A logical scalar. If |
invert |
A logical scalar. If |
removeAllNApos |
A logical scalar. If |
assayNameNA |
A character scalar corresponding to the name of a
read-level assay (if |
A subset RangedSummarizedExperiment
object.
Michael Stadler, Charlotte Soneson
library(SummarizedExperiment) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(modbamfiles, "chr1:6940000-6955000", "a", BPPARAM = BiocParallel::SerialParam()) lapply(assay(se, "mod_prob"), colnames) # subset by read identifiers seSub <- subsetReads(se, c("s1-233e48a7-f379-4dcf-9270-958231125563", "s2-034b625e-6230-4f8d-a713-3a32cd96c298")) lapply(assay(seSub, "mod_prob"), colnames) # subset by a list of indices seSub <- subsetReads(se, list(s1 = c(1, 3), s2 = c(TRUE, FALSE))) lapply(assay(seSub, "mod_prob"), colnames)library(SummarizedExperiment) modbamfiles <- system.file("extdata", c("6mA_1_10reads.bam", "6mA_2_10reads.bam"), package = "SingleMoleculeGenomicsIO") se <- readModBam(modbamfiles, "chr1:6940000-6955000", "a", BPPARAM = BiocParallel::SerialParam()) lapply(assay(se, "mod_prob"), colnames) # subset by read identifiers seSub <- subsetReads(se, c("s1-233e48a7-f379-4dcf-9270-958231125563", "s2-034b625e-6230-4f8d-a713-3a32cd96c298")) lapply(assay(seSub, "mod_prob"), colnames) # subset by a list of indices seSub <- subsetReads(se, list(s1 = c(1, 3), s2 = c(TRUE, FALSE))) lapply(assay(seSub, "mod_prob"), colnames)