Significance level with multiple pairwise comparisons

[SrenBlikdal]

Dear Keegan,

Thank you for a great tool!

I have a question about interpreting the dmrseq results and defining significant regions for downstream analysis when doing multiple pairwise comparisons:

I have compared 5 groups and found DMRs among the 10 pairwise comparisons using dmrseq() and default settings.
When examining the most significant DMRs between the groups eg. the groups D1 and U2 I get:

head(regions_D1_vs_U2) GRanges object with 6 ranges and 7 metadata columns: seqnames ranges strand | L area beta stat pval qval index

| [1] chr03 2145161-2146450 * | 168 27.3299 -0.487466 -34.7782 5.56273e-06 0.0969475 702869-703036 [2] chr12 3569697-3570576 * | 82 27.1275 -0.868408 -27.5666 5.56273e-06 0.0969475 5267243-5267324 **[3] chr19 1148612-1149387 * | 47 33.2763 -1.742097 -27.2877 5.56273e-06 0.0969475 8373117-8373163** [4] chr12 13445474-13446507 * | 173 43.2373 -0.726438 -22.0255 2.22509e-05 0.0969475 5496562-54... [5] chr14 7463232-7463826 * | 61 36.4485 -1.515688 -21.9971 2.22509e-05 0.0969475 6313405-6313465 [6] chr12 4767778-4769403 * | 81 35.3088 1.025720 20.3324 2.22509e-05 0.0969475 5299005-5299085 ------- length(regions_D1_vs_U2) [1] 60630

While the most significant DMRs between D1 and R1 are:

head(regions_D1_vs_R1,20) GRanges object with 20 ranges and 7 metadata columns: seqnames ranges strand | L area beta stat pval qval index

| [1] chr12 3569574-3570615 * | 85 30.8903 -0.993484 -32.8400 5.67202e-06 0.0221645 5178086-5178170 [2] chr21 1851949-1855576 * | 200 104.6914 -1.282397 -29.7122 5.67202e-06 0.0221645 9002818-9003017 [3] chr12 4767624-4770219 * | 104 49.8623 1.134925 26.4022 5.67202e-06 0.0221645 5209353-5209456 [4] chr09 3308515-3309414 * | 55 34.8944 -1.649020 -26.1033 5.67202e-06 0.0221645 3742501-3742555 [5] chr17 14097041-14097874 * | 56 25.2904 -1.145939 -25.1792 5.67202e-06 0.0221645 7784848-7784.. ... ... ... ... . ... ... ... ... ... ... ... [16] chr05 16970459-16971693 * | 56 21.3600 -0.920008 -18.3395 5.67202e-06 0.0221645 2071338-2071.. **[17] chr21 15386084-15386637 * | 25 9.2377 -1.012059 -18.1877 1.13440e-05 0.0262690 9339935-9339..** [18] chr11 3799201-3800276 * | 37 13.5825 -1.040250 -18.1378 1.13440e-05 0.0262690 4748892-4748928 [19] chr04 1652215-1652816 * | 67 23.1132 -0.945565 -18.0487 1.13440e-05 0.0262690 1266843-1266909 [20] chr07 5204199-5207006 * | 89 27.7387 -0.742104 -17.9870 1.13440e-05 0.0262690 2870792-2870880 ------- length(regions_D1_vs_R1) [1] 62523

If I filter all my results for DMRs with a qval < 0.05, I will find that region [3] in D1_vs_U2 is insignificant while eg. region [17] in D1_vs_R1 is significant. To me it seems paradoxical as the pval is lower-, the absolute stat is higher-, the beta difference is higher- and L is higher for region [3] in D1_vs_U2 compared to region [17] in D1_vs_R1.

1) Is my interpretation of a greater difference between region [3] in D1_vs_U2 than region [17] in D1_vs_R1 correct, or did I miss something?

2) Is there any approach that would allow me to filter using a universal significance level and enable me to compare eg the number of DMRs among multiple pairwise comparisons?

Thank you in advance! Best, Soren

[kdkorthauer]

Hi Soren,

Thanks for reaching out. I'm happy to give a few pointers to help interpret these results. Please see my responses below.

1. Is my interpretation of a greater difference between region [3] in D1_vs_U2 than region [17] in D1_vs_R1 correct, or did I miss something?

You are correct that the estimate of the effect size for the difference between D1 and U2 in region 3 is greater than the effect size estimate for the difference between D1 and R1 in region 17. In addition, the marginal p-value for the first is smaller than the second. However, it is important to note that you are applying 10 separate pairwise comparisons. In each one, a comparison-specific null distribution is generated. In addition, the adjustment for multiple comparisons is done separately in each comparison, which each has a different number of candidate regions, and signal to noise ratio. So this means that each comparison has its own mapping from raw p-values to FDR estimates (q-values). This means that comparing raw p-values across comparisons after computing FDR values for each comparison separately is not meaningful.

2. Is there any approach that would allow me to filter using a universal significance level and enable me to compare eg the number of DMRs among multiple pairwise comparisons?

One approach you could take is to take the entire collection of raw p-values from all candidate regions from all pairwise comparisons, and compute FDR estimates (e.g. using the p.adjust function). Another approach you may consider if you are equally interested in differences among any of the 5 groups is to try analyzing the entire dataset with one run of dmrseq, where instead of the covariate of interest being a two group factor, it is a 5 level factor. This is akin to an ANOVA like test where you are testing the null hypothesis that all groups are equal vs the alternative that at least one group is different than the rest. However, if you are interested in drilling down to pinpoint exactly which groups the differences were present, then the pairwise tests would still be necessary.

Hope that helps!

Best, Keegan

[SrenBlikdal]

Hi Keegan,

Thank you very much for this quick and elaborate answer, a great help indeed! Using all the raw p-values to compute FDR estimates, definitely seems like a great way of proceeding.

I will close this issue.

Best, Soren