ci: unit test for --no_collapse, drop irrelevant comment

docs/content/tools/alignmentSieve.rst

@@ -78,9 +78,9 @@ and the ``-1 4`` set would produce the following::
 
          ------------------------ fragment
 
-             --------------------- shifted fragment
+    -------------------------------- shifted fragment
 
-As can be seen, such fragments are considered to be on the ``-`` strand, so negative values then shift to the left on its frame of reference (thus, to the right relative to the ``+`` strand).
+As can be seen, such fragments are considered to be on the ``-`` strand, so negative values then shift to the left on its frame of reference (thus, to the right relative to the ``+`` strand). Note that, just like the ``-5 3`` example above, both ends of the fragment move outward (the left end shifts left by 4, the right end shifts right by 1): the sign of each value only determines its direction, not whether the two ends move toward or away from each other.
 
 .. note::
     If the ``--shift`` or ``--ATACshift`` options are used, then only properly-paired reads will be used.

docs/content/tools/bamCompare.rst

@@ -29,3 +29,5 @@ The basic algorithm works proceeds in two steps:
    :nodefault:
 
 .. note:: As of deepTools 4.0.0, ``bamCompare`` uses a new Rust-backed core. ``--blackListFileName`` may be gzip-compressed and blacklist filtering is done at base-pair resolution rather than by rejecting whole genomic chunks. The SES scaling method and ``--ignoreDuplicates`` have both been removed; for duplicate removal use ``--samFlagExclude`` against a BAM file with duplicates marked.
+
+.. note:: ``--normalizeUsing RPGC`` is **not supported** by ``bamCompare`` and will exit with an error. Use ``--normalizeUsing RPKM``, ``CPM`` or ``BPM``, or leave ``--scaleFactorsMethod`` at its default (``readCount``) to equalize sequencing depth between the two samples instead.

pydeeptools/deeptools/alignmentSieve2.py

@@ -79,8 +79,12 @@ def parseArguments():
     filtering = parser.add_argument_group('Optional arguments')
 
     filtering.add_argument('--filterRNAstrand',
-                           help='Selects RNA-seq reads (single-end or paired-end) in '
-                                'the given strand. (Default: %(default)s)',
+                           help='Selects RNA-seq reads (single-end or paired-end) originating from genes '
+                                'on the given strand. This option assumes a standard dUTP-based library '
+                                'preparation (that is, --filterRNAstrand=forward keeps minus-strand reads, '
+                                'which originally came from genes on the forward strand using a dUTP-based '
+                                'method). Consider using --samFlagExclude instead for filtering by strand in '
+                                'other contexts. (Default: %(default)s)',
                            choices=['forward', 'reverse', 'None'],
                            default='None')
 

pydeeptools/deeptools/bamCompare2.py

@@ -150,6 +150,11 @@ def getOptionalArgs():
 def process_args(args=None):
     args = parseArguments().parse_args(args)
 
+    if args.normalizeUsing == "RPGC":
+        sys.exit(
+            "RPGC normalization (--normalizeUsing RPGC) is not supported with bamCompare. "
+        )
+
     if args.smoothLength and args.smoothLength <= args.binSize:
         print(f"Warning: the smooth length given ({args.smoothLength}) is smaller than the bin "
               f"size ({args.binSize}).\n\n No smoothing will be "

pydeeptools/deeptools/correlation.py

@@ -38,10 +38,12 @@ def __init__(self, matrix_file,
         self.corr_matrix = None  # correlation matrix
         self.column_order = None
         if labels is not None:
-            # test that the length of labels
-            # corresponds to the length of
-            # samples
-
+            if len(labels) != self.matrix.shape[1]:
+                sys.exit(
+                    f"\nThe number of labels provided ({len(labels)}) does not match the number "
+                    f"of samples in the matrix ({self.matrix.shape[1]}). Please provide exactly one label "
+                    "per sample.\n"
+                )
             self.labels = labels
         self.labels = [toString(x) for x in self.labels]
 

pydeeptools/deeptools/countReadsPerBin.py

@@ -898,7 +898,8 @@ def get_fragment_from_read(self, read):
                                             f"end for read {read.query_name}"
         return [(fragmentStart, fragmentEnd)]
 
-    def getSmoothRange(self, tileIndex, tileSize, smoothRange, maxPosition):
+    @staticmethod
+    def getSmoothRange(tileIndex, tileSize, smoothRange, maxPosition):
         """
         Given a tile index position and a tile size (length), return the a new indices
         over a larger range, called the smoothRange.

pydeeptools/deeptools/estimateReadFiltering.py

@@ -99,8 +99,12 @@ def parseArguments():
     filtering = parser.add_argument_group('Optional arguments')
 
     filtering.add_argument('--filterRNAstrand',
-                           help='Selects RNA-seq reads (single-end or paired-end) in '
-                                'the given strand. (Default: %(default)s)',
+                           help='Selects RNA-seq reads (single-end or paired-end) originating from genes '
+                                'on the given strand. This option assumes a standard dUTP-based library '
+                                'preparation (that is, --filterRNAstrand=forward keeps minus-strand reads, '
+                                'which originally came from genes on the forward strand using a dUTP-based '
+                                'method). Consider using --samFlagExclude instead for filtering by strand in '
+                                'other contexts. (Default: %(default)s)',
                            choices=['forward', 'reverse'],
                            default=None)
 

pydeeptools/deeptools/parserCommon.py

@@ -231,7 +231,8 @@ def normalization_options():
 
     group.add_argument('--normalizeUsing',
                        help='Use one of the entered methods to '
-                       'normalize the number of reads per bin. By default, no normalization is performed. '
+                       'normalize the number of reads per bin.'
+                       'None = the default and equivalent to not setting this option at all. '
                        'RPKM = Reads Per Kilobase per Million mapped reads; '
                        'CPM = Counts Per Million mapped reads, same as CPM in RNA-seq; '
                        'BPM = Bins Per Million mapped reads, same as TPM in RNA-seq; '
@@ -245,7 +246,6 @@ def normalization_options():
                        'sum of all reads per bin (in millions). '
                        'RPGC (per bin) = number of reads per bin / '
                        'scaling factor for 1x average coverage. '
-                       'None = the default and equivalent to not setting this option at all. '
                        'This scaling factor, in turn, is determined from the '
                        'sequencing depth: (total number of mapped reads * fragment length) / '
                        'effective genome size.\nThe scaling factor used '

pydeeptools/deeptools/test/test_bamCoverage_and_bamCompare.py

@@ -2,6 +2,8 @@
 import tempfile
 from os import unlink
 
+import pytest
+
 import deeptools.bamCompare2 as bam_comp
 import deeptools.bamCoverage2 as bam_cov
 
@@ -126,22 +128,6 @@ def test_bam_coverage_scaleFactor():
         assert f"{resp}" == f"{expected}", f"{resp} != {expected}"
         unlink(outfile)
 
-
-# def test_bam_coverage_filtering():
-#     _, outfile = tempfile.mkstemp(suffix=".bg")
-#     #for fname in [BAMFILE_B, CRAMFILE_B]:
-#     for fname in [BAMFILE_B]:
-#         args = "--bam {} -o {} --outFileFormat bedgraph --ignoreDuplicates --verbose".format(fname, outfile).split()
-#         bam_cov.main(args)
-
-#         _foo = open(outfile, 'r')
-#         resp = _foo.readlines()
-#         _foo.close()
-#         expected = ['3R\t0\t50\t0\n', '3R\t50\t200\t1\n']
-#         assert resp == expected, "{} != {}".format(resp, expected)
-#         unlink(outfile)
-
-
 def test_bam_compare_arguments():
     """
     Test minimal command line args for bamCoverage. The ratio
@@ -162,6 +148,16 @@ def test_bam_compare_arguments():
         unlink(outfile)
 
 
+def test_bam_compare_rpgc_rejected():
+    """
+    RPGC normalization in bamcompare -> failure.
+    """
+    args = f"--bamfile1 {BAMFILE_A} --bamfile2 {BAMFILE_B} --normalizeUsing RPGC " \
+           "--effectiveGenomeSize 1000000 -o /dev/null".split()
+    with pytest.raises(SystemExit):
+        bam_comp.main(args)
+
+
 def test_bam_compare_diff_files():
     """
     Test with two different files
@@ -212,38 +208,6 @@ def test_bam_compare_ZoverZ():
     unlink(outfile)
 
 
-# def test_get_num_kept_reads():
-#     """
-#     Test the scale factor functions
-#     """
-#     for fname in [BAMFILE_A, CRAMFILE_A]:
-#         args = "--bam {}  -o /tmp/test".format(fname).split()
-
-#         args = bam_cov.process_args(args)
-#         num_kept_reads, total_reads = gs.get_num_kept_reads(args, None)
-
-#         # bam file 1 has 2 reads in 3R and 2 read in chr_cigar
-#         assert num_kept_reads == 3, "num_kept_reads is wrong"
-#         assert total_reads == 3, "num total reads is wrong"
-
-#         # ignore chr_cigar to count the total number of reads
-#         args = "--bam {} --ignoreForNormalization chr_cigar  -o /tmp/test".format(fname).split()
-#         args = bam_cov.process_args(args)
-#         num_kept_reads, total_reads = gs.get_num_kept_reads(args, None)
-
-#         # the  number of kept reads should be 2 as the read on chr_cigar is skipped
-#         assert num_kept_reads == 2, "num_kept_reads is wrong ({})".format(num_kept_reads)
-
-#         # test filtering by read direction. Only forward reads are kept
-#         args = "--bam {}  -o /tmp/test --samFlagExclude 16 --ignoreForNormalization chr_cigar ".format(fname).split()
-
-#         args = bam_cov.process_args(args)
-#         num_kept_reads, total_reads = gs.get_num_kept_reads(args, None)
-
-#         # only one forward read is expected in
-#         assert num_kept_reads == 1, "num_kept_reads is wrong"
-
-
 def test_bam_compare_diff_files_skipnas():
     """
     Test skipnas
@@ -585,3 +549,67 @@ def test_bam_compare_filter_blacklist():
         ]
         assert f"{resp}" == f"{expected}", f"{resp} != {expected}"
         unlink(outfile)
+
+def test_bam_coverage_nocollapse():
+    """
+    Test --no_collapse in bamcoverage
+    """
+    _, outfile = tempfile.mkstemp(suffix=".bg")
+    args = f"-b {BAMFILE_A} -of bedgraph -bs 20 --no_collapse -o {outfile}"
+    print(args)
+    args = args.split()
+    bam_cov.main(args)
+
+    with open(outfile, 'r') as _foo:
+        resp = _foo.readlines()
+        expected = [
+            "3R\t0\t20\t0\n",
+            "3R\t20\t40\t0\n",
+            "3R\t40\t60\t0\n",
+            "3R\t60\t80\t0\n",
+            "3R\t80\t100\t0\n",
+            "3R\t100\t120\t1\n",
+            "3R\t120\t140\t1\n",
+            "3R\t140\t160\t2\n",
+            "3R\t160\t180\t1\n",
+            "3R\t180\t200\t1\n",
+            "chr_cigar\t0\t20\t1\n",
+            "chr_cigar\t20\t40\t1\n",
+            "chr_cigar\t40\t60\t1\n",
+            "chr_cigar\t60\t80\t0\n",
+            "chr_cigar\t80\t100\t0\n",
+            "chr_cigar\t100\t120\t0\n",
+            "chr_cigar\t120\t140\t0\n",
+            "chr_cigar\t140\t160\t0\n",
+            "chr_cigar\t160\t180\t0\n",
+            "chr_cigar\t180\t200\t0\n",
+        ]
+        assert f"{resp}" == f"{expected}", f"{resp} != {expected}"
+        unlink(outfile)
+
+def test_bam_compare_nocollapse():
+    """
+    Test --no_collapse in bamcompare
+    """
+    _, outfile = tempfile.mkstemp(suffix=".bg")
+    args = f"-b1 {BAMFILE_A} -b2 {BAMFILE_B} -of bedgraph -bs 20 --no_collapse -o {outfile}"
+    print(args)
+    args = args.split()
+    bam_comp.main(args)
+
+    with open(outfile, 'r') as _foo:
+        resp = _foo.readlines()
+        expected = [
+            "3R\t0\t20\t0\n",
+            "3R\t20\t40\t0\n",
+            "3R\t40\t60\t-0.58\n",
+            "3R\t60\t80\t-0.58\n",
+            "3R\t80\t100\t-0.58\n",
+            "3R\t100\t120\t0.42\n",
+            "3R\t120\t140\t0.42\n",
+            "3R\t140\t160\t0.26\n",
+            "3R\t160\t180\t0\n",
+            "3R\t180\t200\t0\n",
+        ]
+        assert f"{resp}" == f"{expected}", f"{resp} != {expected}"
+        unlink(outfile)

pydeeptools/deeptools/test/test_plotcorrelation.py

@@ -15,6 +15,14 @@
 COR_PLOT_GG_2 = ROOT + "plotCorrelation_result2_ggplot.png"
 
 
+def test_correlation_labels_length_mismatch_exits():
+    _, out_png = tempfile.mkstemp(suffix=".png")
+    args = f"--corData {COR_DATA_IN1} -p heatmap -c pearson -o {out_png} " \
+           "--labels sample1 sample2 sample3".split()
+    with pytest.raises(SystemExit):
+        pc.main(args)
+
+
 @pytest.mark.filterwarnings(
     "ignore:Attempting to set identical low and high xlims:UserWarning"
 )

pydeeptools/deeptools/writeBedGraph_bam_and_bw.py

@@ -10,6 +10,7 @@
 
 # own module
 from deeptools import bamHandler, mapReduce
+from deeptools.countReadsPerBin import CountReadsPerBin
 from deeptools.utilities import getCommonChrNames, toBytes
 from deeptools.writeBedGraph import *
 
@@ -81,7 +82,7 @@ def writeBedGraph_worker(
             tileCoverage = []
             for index in range(len(bamOrBwFileList)):
                 if smoothLength > 0:
-                    vectorStart, vectorEnd = getSmoothRange(
+                    vectorStart, vectorEnd = CountReadsPerBin.getSmoothRange(
                         tileIndex, tileSize, smoothLength, lengthCoverage)
                     tileCoverage.append(
                         np.mean(coverage[index][vectorStart:vectorEnd]))