Fit GLM and return DMR data in one step

Calculate one or more contrasts to find Differentially Methylated Regions (DMRs).

calculateDMRs(
  qseaSet,
  variable = NULL,
  covariates = NULL,
  contrasts = NULL,
  minReadCount = 0,
  minNRPM = 1,
  checkPVals = TRUE,
  FDRthres = 0.05,
  keepPvals = FALSE,
  formula = NULL,
  keepContrastMeans = TRUE,
  keepData = FALSE,
  keepGroupMeans = FALSE,
  direction = "both",
  calcDispersionAll = FALSE
)

Arguments

qseaSet

qseaSet. Input object containing counts and sample metadata.

variable

character(1) or NULL. Sample-table column used as the primary explanatory variable for contrasts. Default: NULL (must be supplied).

covariates

character() or NULL. Additional covariate column names (e.g., batch, patient). Default: NULL.

contrasts

data.frame, character(1), or NULL. Contrast specification. If a data frame, must contain columns group1 and group2. If a string, supports:

"A_vs_B" — a single explicit contrast,
"All" / "all" — all pairwise contrasts among levels of variable,
"All_vs_X" — each level vs X,
"X_vs_All" — X vs each other level,
"first" — only the first possible contrast (based on factor order). Default: NULL (fit the GLM without adding contrasts).

minReadCount

numeric(1). Minimum raw count required in at least one sample to keep a window. Default: 0.

minNRPM

numeric(1). Minimum NRPM required in at least one sample to keep a window. Default: 1.

checkPVals

logical(1). If TRUE, stop/warn when an excessive proportion of p-values are exactly zero (instability check). Default: TRUE.

FDRthres

numeric(1). False discovery rate threshold for calling significance. Default: 0.05.

keepPvals

logical(1). If TRUE, include raw (unadjusted) p-values in the output. Default: FALSE.

formula

formula or NULL. Optional model formula overriding variable/covariates. Default: NULL.

keepContrastMeans

logical(1). If FALSE, drop contrast mean columns from the output. Default: TRUE.

keepData

logical(1). If TRUE, include per-sample data columns for significant windows. Default: FALSE.

keepGroupMeans

logical(1). If TRUE, include group mean columns (based on variable). Default: FALSE.

direction

character(1). Direction for calling significance, passed to qsea::isSignificant(): one of "up", "down", or "both". Default: "both".

calcDispersionAll

logical(1). If TRUE, samples not present in any contrast still contribute to the initial GLM fit/dispersion estimation (so adding samples can change DMRs). Default: FALSE.

Value

A tibble with one row per significant region/window containing:

per-contrast statistics (e.g., log2FC, deltaBeta, test statistic),
multiple-testing results (FDR; optionally raw p-values if keepPvals),
optional group means (if keepGroupMeans = TRUE),
optional per-sample columns (if keepData = TRUE),
metadata columns for genomic coordinates and window attributes.

Details

Contrast strings are parsed as follows:

"A_vs_B" creates a single contrast A vs B;
"All"/"all" expands to all pairwise contrasts among levels of variable (via makeAllContrasts());
"All_vs_X" creates each level vs X;
"X_vs_All" creates X vs each other level;
"first" constructs only the first possible pair (based on factor order). Internally, this function:

filters windows by minReadCount/minNRPM, 2) fits a Negative Binomial Generalised Linear Model including dispersion estimates, 3) adds contrasts as requested, calculating adjusted p-values (adjPval), log2 fold changes (log2FC) and differences in beta values (deltaBeta), 4) creates an output table showing any window which is differentially methylated (with adjusted p-value below FDRthres) in at least one contrast. Additional columns will be added containing data values if keepData and keepGroupMeans are true.

Examples

data(exampleTumourNormal, package = "mesa")

# A single explicit contrast
exampleTumourNormal %>%
    calculateDMRs(variable = "type", contrasts = "LUAD_vs_NormalLung")
#> Generating table of nrpm values for 819 regions across 5 samples.
#> Fitting initial GLM on 255 windows, without using parallelisation.
#> Performing contrast LUAD_vs_NormalLung
#> Contrasts Complete
#> selecting regions from contrast LUAD_vs_NormalLung
#> selected 46/255 windows
#> adding test results from LUAD_vs_NormalLung
#> # A tibble: 46 × 17
#>    seqnames    start      end CpG_density LUAD_vs_NormalLung_log2FC
#>    <fct>       <int>    <int>       <dbl>                     <dbl>
#>  1 7        25107301 25107600        8.13                     -1.71
#>  2 7        27096001 27096300       18.8                       3.72
#>  3 7        27096301 27096600       19.9                       3.38
#>  4 7        27096601 27096900       12.6                       2.44
#>  5 7        27096901 27097200       17.9                       2.99
#>  6 7        27098401 27098700       12.1                       1.11
#>  7 7        27103501 27103800       13.1                       1.72
#>  8 7        27106501 27106800       19.0                       1.06
#>  9 7        27106801 27107100       13.1                       1.38
#> 10 7        27108301 27108600       15.7                       1.06
#> # ℹ 36 more rows
#> # ℹ 12 more variables: LUAD_vs_NormalLung_adjPval <dbl>,
#> #   LUAD_vs_NormalLung_deltaBeta <dbl>, NormalColon_beta_means <dbl>,
#> #   CRC_beta_means <dbl>, NormalLung_beta_means <dbl>, LUAD_beta_means <dbl>,
#> #   LUSC_beta_means <dbl>, NormalColon_nrpm_means <dbl>, CRC_nrpm_means <dbl>,
#> #   NormalLung_nrpm_means <dbl>, LUAD_nrpm_means <dbl>, LUSC_nrpm_means <dbl>

# All pairwise contrasts among levels of 'type'
exampleTumourNormal %>%
    calculateDMRs(variable = "type", contrasts = "all")
#> Calculating 10 possible contrasts on the type column.
#> Generating table of nrpm values for 819 regions across 10 samples.
#> Fitting initial GLM on 309 windows, without using parallelisation.
#> Performing contrast CRC_vs_LUAD
#> Performing contrast CRC_vs_LUSC
#> Performing contrast CRC_vs_NormalColon
#> Performing contrast CRC_vs_NormalLung
#> Performing contrast LUAD_vs_LUSC
#> Performing contrast LUAD_vs_NormalColon
#> Performing contrast LUAD_vs_NormalLung
#> Performing contrast LUSC_vs_NormalColon
#> Performing contrast LUSC_vs_NormalLung
#> Performing contrast NormalColon_vs_NormalLung
#> Contrasts Complete
#> selecting regions from contrast CRC_vs_LUAD
#> selected 33/309 windows
#> selecting regions from contrast CRC_vs_LUSC
#> selected 38/309 windows
#> selecting regions from contrast CRC_vs_NormalColon
#> selected 1/309 windows
#> selecting regions from contrast CRC_vs_NormalLung
#> selected 44/309 windows
#> selecting regions from contrast LUAD_vs_LUSC
#> selected 12/309 windows
#> selecting regions from contrast LUAD_vs_NormalColon
#> selected 13/309 windows
#> selecting regions from contrast LUAD_vs_NormalLung
#> selected 23/309 windows
#> selecting regions from contrast LUSC_vs_NormalColon
#> selected 38/309 windows
#> selecting regions from contrast LUSC_vs_NormalLung
#> selected 57/309 windows
#> selecting regions from contrast NormalColon_vs_NormalLung
#> selected 18/309 windows
#> adding test results from CRC_vs_LUAD
#> adding test results from CRC_vs_LUSC
#> adding test results from CRC_vs_NormalColon
#> adding test results from CRC_vs_NormalLung
#> adding test results from LUAD_vs_LUSC
#> adding test results from LUAD_vs_NormalColon
#> adding test results from LUAD_vs_NormalLung
#> adding test results from LUSC_vs_NormalColon
#> adding test results from LUSC_vs_NormalLung
#> adding test results from NormalColon_vs_NormalLung
#> # A tibble: 109 × 44
#>    seqnames    start      end CpG_density CRC_vs_LUAD_log2FC CRC_vs_LUAD_adjPval
#>    <fct>       <int>    <int>       <dbl>              <dbl>               <dbl>
#>  1 7        25018801 25019100        15.3              -1.99          0.0496    
#>  2 7        25852801 25853100        18.4               2.50          0.000410  
#>  3 7        25860601 25860900        11.6               2.35          0.0115    
#>  4 7        27087901 27088200        18.7               1.72          0.00916   
#>  5 7        27088201 27088500        19.7               2.01          0.000575  
#>  6 7        27096001 27096300        18.8              -2.74          0.00250   
#>  7 7        27096301 27096600        19.9              -4.15          0.00000512
#>  8 7        27096601 27096900        12.6              -2.97          0.0244    
#>  9 7        27096901 27097200        17.9              -4.04          0.00806   
#> 10 7        27102001 27102300        17.1               1.82          0.00392   
#> # ℹ 99 more rows
#> # ℹ 38 more variables: CRC_vs_LUAD_deltaBeta <dbl>, CRC_vs_LUSC_log2FC <dbl>,
#> #   CRC_vs_LUSC_adjPval <dbl>, CRC_vs_LUSC_deltaBeta <dbl>,
#> #   CRC_vs_NormalColon_log2FC <dbl>, CRC_vs_NormalColon_adjPval <dbl>,
#> #   CRC_vs_NormalColon_deltaBeta <dbl>, CRC_vs_NormalLung_log2FC <dbl>,
#> #   CRC_vs_NormalLung_adjPval <dbl>, CRC_vs_NormalLung_deltaBeta <dbl>,
#> #   LUAD_vs_LUSC_log2FC <dbl>, LUAD_vs_LUSC_adjPval <dbl>, …

# Tighter FDR threshold
exampleTumourNormal %>%
    calculateDMRs(variable = "type", contrasts = "all", FDRthres = 1e-4)
#> Calculating 10 possible contrasts on the type column.
#> Generating table of nrpm values for 819 regions across 10 samples.
#> Fitting initial GLM on 309 windows, without using parallelisation.
#> Performing contrast CRC_vs_LUAD
#> Performing contrast CRC_vs_LUSC
#> Performing contrast CRC_vs_NormalColon
#> Performing contrast CRC_vs_NormalLung
#> Performing contrast LUAD_vs_LUSC
#> Performing contrast LUAD_vs_NormalColon
#> Performing contrast LUAD_vs_NormalLung
#> Performing contrast LUSC_vs_NormalColon
#> Performing contrast LUSC_vs_NormalLung
#> Performing contrast NormalColon_vs_NormalLung
#> Contrasts Complete
#> selecting regions from contrast CRC_vs_LUAD
#> selected 2/309 windows
#> selecting regions from contrast CRC_vs_LUSC
#> selected 7/309 windows
#> selecting regions from contrast CRC_vs_NormalColon
#> selected 0/309 windows
#> selecting regions from contrast CRC_vs_NormalLung
#> selected 9/309 windows
#> selecting regions from contrast LUAD_vs_LUSC
#> selected 0/309 windows
#> selecting regions from contrast LUAD_vs_NormalColon
#> selected 0/309 windows
#> selecting regions from contrast LUAD_vs_NormalLung
#> selected 2/309 windows
#> selecting regions from contrast LUSC_vs_NormalColon
#> selected 3/309 windows
#> selecting regions from contrast LUSC_vs_NormalLung
#> selected 12/309 windows
#> selecting regions from contrast NormalColon_vs_NormalLung
#> selected 2/309 windows
#> adding test results from CRC_vs_LUAD
#> adding test results from CRC_vs_LUSC
#> adding test results from CRC_vs_NormalColon
#> adding test results from CRC_vs_NormalLung
#> adding test results from LUAD_vs_LUSC
#> adding test results from LUAD_vs_NormalColon
#> adding test results from LUAD_vs_NormalLung
#> adding test results from LUSC_vs_NormalColon
#> adding test results from LUSC_vs_NormalLung
#> adding test results from NormalColon_vs_NormalLung
#> # A tibble: 24 × 44
#>    seqnames    start      end CpG_density CRC_vs_LUAD_log2FC CRC_vs_LUAD_adjPval
#>    <fct>       <int>    <int>       <dbl>              <dbl>               <dbl>
#>  1 7        27096301 27096600       19.9              -4.15           0.00000512
#>  2 7        27130201 27130500       28.2              -1.88           0.0000619 
#>  3 7        27151201 27151500       16.5              -3.68           0.00292   
#>  4 7        27151501 27151800        8.60             -1.26           0.587     
#>  5 7        27152401 27152700       15.6              -1.51           0.328     
#>  6 7        27165601 27165900       13.9              -3.10           0.00953   
#>  7 7        27166201 27166500       20.0              -3.15           0.00275   
#>  8 7        27168601 27168900        8.79             -1.66           0.326     
#>  9 7        27185401 27185700       13.6              -0.972          0.657     
#> 10 7        25860601 25860900       11.6               2.35           0.0115    
#> # ℹ 14 more rows
#> # ℹ 38 more variables: CRC_vs_LUAD_deltaBeta <dbl>, CRC_vs_LUSC_log2FC <dbl>,
#> #   CRC_vs_LUSC_adjPval <dbl>, CRC_vs_LUSC_deltaBeta <dbl>,
#> #   CRC_vs_NormalColon_log2FC <dbl>, CRC_vs_NormalColon_adjPval <dbl>,
#> #   CRC_vs_NormalColon_deltaBeta <dbl>, CRC_vs_NormalLung_log2FC <dbl>,
#> #   CRC_vs_NormalLung_adjPval <dbl>, CRC_vs_NormalLung_deltaBeta <dbl>,
#> #   LUAD_vs_LUSC_log2FC <dbl>, LUAD_vs_LUSC_adjPval <dbl>, …

# Return all per-sample columns for the first available contrast on 'tumour'
exampleTumourNormal %>%
    calculateDMRs(variable = "tumour", contrasts = "first", keepData = TRUE)
#> Calculating the first possible contrast on the tumour column.
#> Generating table of nrpm values for 819 regions across 10 samples.
#> Fitting initial GLM on 309 windows, without using parallelisation.
#> Performing contrast Normal_vs_Tumour
#> Contrasts Complete
#> selecting regions from contrast Normal_vs_Tumour
#> selected 23/309 windows
#> adding test results from Normal_vs_Tumour
#> # A tibble: 23 × 31
#>    seqnames    start      end CpG_density Normal_vs_Tumour_log2FC
#>    <fct>       <int>    <int>       <dbl>                   <dbl>
#>  1 7        27102901 27103200       10.8                   -1.94 
#>  2 7        27103501 27103800       13.1                   -1.88 
#>  3 7        27103801 27104100        7.83                  -1.33 
#>  4 7        27108301 27108600       15.7                   -1.18 
#>  5 7        27108601 27108900       16.1                   -0.933
#>  6 7        27144601 27144900       11.1                   -1.47 
#>  7 7        27144901 27145200       16.3                   -1.30 
#>  8 7        27156301 27156600       18.5                   -2.57 
#>  9 7        27159901 27160200        9.48                   1.16 
#> 10 7        27164701 27165000       13.8                   -2.25 
#> # ℹ 13 more rows
#> # ℹ 26 more variables: Normal_vs_Tumour_adjPval <dbl>,
#> #   Normal_vs_Tumour_deltaBeta <dbl>, Colon1_N_beta <dbl>, Colon1_T_beta <dbl>,
#> #   Colon2_N_beta <dbl>, Colon2_T_beta <dbl>, Lung1_N_beta <dbl>,
#> #   Lung1_T_beta <dbl>, Lung2_N_beta <dbl>, Lung2_T_beta <dbl>,
#> #   Lung3_N_beta <dbl>, Lung3_T_beta <dbl>, Colon1_N_nrpm <dbl>,
#> #   Colon1_T_nrpm <dbl>, Colon2_N_nrpm <dbl>, Colon2_T_nrpm <dbl>, …

Fit GLM and return DMR data in one step

Arguments

Value

Details

See also

Examples