Adding the possibility to use a threshold if no multitesting correcti…

…on is applied. See issue #145
deeptools · Nov 8, 2017 · 88be5eb · 88be5eb
1 parent a2b5f81
commit 88be5eb
Show file tree

Hide file tree

Showing 5 changed files with 334 additions and 334 deletions.
diff --git a/hicexplorer/hicFindTADs.py b/hicexplorer/hicFindTADs.py
@@ -131,7 +131,9 @@ def parse_arguments(args=None):
                              'is estimated by comparing the distribution (Wilcoxon ranksum) of '
                              'the  zscores between the left and right '
                              'regions (diamond) at the local minimum with the matrix zscores for a '
-                             'diamond at --minDepth to the left and a diamond --minDepth to the right. ',
+                             'diamond at --minDepth to the left and a diamond --minDepth to the right. '
+                             'If --correctForMultipleTesting is \'None\' the threshold is applied on the '
+                             'raw p-values without any multiple testing correction. Set it to \'1\' if no threshold should be used.',
                         type=float,
                         default=0.01)
 
@@ -934,18 +936,19 @@ def save_domains_and_boundaries(self, prefix):
             elif self.correct_for_multiple_testing == 'bonferroni':
                 if delta_of_min[idx] >= self.delta and idx in pvalue_of_min and pvalue_of_min[idx] <= self.threshold_comparisons:
                     filtered_min_idx += [idx]
-            else:
-                if delta_of_min[idx] >= self.delta:
+            elif self.correct_for_multiple_testing == 'None':
+                if delta_of_min[idx] >= self.delta and idx in pvalue_of_min and pvalue_of_min[idx] <= self.threshold_comparisons:
                     filtered_min_idx += [idx]
 
         if self.correct_for_multiple_testing == 'fdr':
-            log.info("Number of boundaries for delta {}, qval {}: {}".format(self.delta, self.threshold_comparisons,
+            log.info("FDR correction. Number of boundaries for delta {}, qval {}: {}".format(self.delta, self.threshold_comparisons,
                                                                              len(filtered_min_idx)))
         elif self.correct_for_multiple_testing == 'bonferroni':
-            log.info("Number of boundaries for delta {} and pval {}: {}".format(self.delta, self.threshold_comparisons,
+            log.info("Bonferroni correction. Number of boundaries for delta {} and pval {}: {}".format(self.delta, self.threshold_comparisons,
                                                                                 len(filtered_min_idx)))
         else:
-            log.info("Number of boundaries for delta {}: {}".format(self.delta, len(filtered_min_idx)))
+            log.info("No multiple testing correction. Number of boundaries for delta {}: {}, used threshold: {}".format(self.delta, len(filtered_min_idx), self.threshold_comparisons))
+
         count = 1
         with open(prefix + '_boundaries.bed', 'w') as file_boundary_bin, open(prefix + '_domains.bed', 'w') as file_domains, open(prefix + '_boundaries.gff', 'w') as gff:
             for idx, min_bin_id in enumerate(filtered_min_idx):

diff --git a/hicexplorer/test/test_data/find_TADs/None/multiNone_boundaries.bed b/hicexplorer/test/test_data/find_TADs/None/multiNone_boundaries.bed
@@ -1080,7 +1080,6 @@ chr2R	20812500	20817500	B18174	-0.074736500000	.
 chr2R	20937500	20942500	B18199	-0.034997750000	.
 chr2R	20972500	20977500	B18206	-0.041723500000	.
 chr2R	21072500	21077500	B18226	0.089617250000	.
-chr3R	47500	52500	B18251	-0.087740250000	.
 chr3R	162500	167500	B18274	-0.033048250000	.
 chr3R	227500	232500	B18287	-0.092860000000	.
 chr3R	292500	297500	B18300	-0.042622250000	.

diff --git a/hicexplorer/test/test_data/find_TADs/None/multiNone_boundaries.gff b/hicexplorer/test/test_data/find_TADs/None/multiNone_boundaries.gff
@@ -1080,7 +1080,6 @@ chr2R	HiCExplorer	boundary	20812500	20817500	-0.074736500000	.	.	ID=B18174;delta
 chr2R	HiCExplorer	boundary	20937500	20942500	-0.034997750000	.	.	ID=B18199;delta=0.018418513158;pvalue=0.876302477922;tad_sep=-0.034997750000
 chr2R	HiCExplorer	boundary	20972500	20977500	-0.041723500000	.	.	ID=B18206;delta=0.048636776316;pvalue=0.654246312620;tad_sep=-0.041723500000
 chr2R	HiCExplorer	boundary	21072500	21077500	0.089617250000	.	.	ID=B18226;delta=0.129759171053;pvalue=0.005235661183;tad_sep=0.089617250000
-chr3R	HiCExplorer	boundary	47500	52500	-0.087740250000	.	.	ID=B18251;delta=0.058089026316;pvalue=nan;tad_sep=-0.087740250000
 chr3R	HiCExplorer	boundary	162500	167500	-0.033048250000	.	.	ID=B18274;delta=0.024394500000;pvalue=0.917163037297;tad_sep=-0.033048250000
 chr3R	HiCExplorer	boundary	227500	232500	-0.092860000000	.	.	ID=B18287;delta=0.049449078947;pvalue=0.416240922379;tad_sep=-0.092860000000
 chr3R	HiCExplorer	boundary	292500	297500	-0.042622250000	.	.	ID=B18300;delta=0.058786013158;pvalue=0.173195074149;tad_sep=-0.042622250000