diff --git a/docs/content/tools/alignmentSieve.rst b/docs/content/tools/alignmentSieve.rst
index 5da71de7ee..f10b165fbd 100644
--- a/docs/content/tools/alignmentSieve.rst
+++ b/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.
diff --git a/docs/content/tools/bamCompare.rst b/docs/content/tools/bamCompare.rst
index 7b673d1272..2e4deae470 100644
--- a/docs/content/tools/bamCompare.rst
+++ b/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.
diff --git a/pydeeptools/deeptools/alignmentSieve2.py b/pydeeptools/deeptools/alignmentSieve2.py
index 44e47bf0d9..39b315b93f 100644
--- a/pydeeptools/deeptools/alignmentSieve2.py
+++ b/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')
 
diff --git a/pydeeptools/deeptools/bamCompare2.py b/pydeeptools/deeptools/bamCompare2.py
index dfc909112f..fceb36822e 100644
--- a/pydeeptools/deeptools/bamCompare2.py
+++ b/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 "
diff --git a/pydeeptools/deeptools/correlation.py b/pydeeptools/deeptools/correlation.py
index 1e59963f40..bcdbc13308 100644
--- a/pydeeptools/deeptools/correlation.py
+++ b/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]
 
diff --git a/pydeeptools/deeptools/countReadsPerBin.py b/pydeeptools/deeptools/countReadsPerBin.py
index 4668befb27..0e7b6f216c 100644
--- a/pydeeptools/deeptools/countReadsPerBin.py
+++ b/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.
diff --git a/pydeeptools/deeptools/estimateReadFiltering.py b/pydeeptools/deeptools/estimateReadFiltering.py
index d44aa5d42d..469e9c10b4 100644
--- a/pydeeptools/deeptools/estimateReadFiltering.py
+++ b/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)
 
diff --git a/pydeeptools/deeptools/parserCommon.py b/pydeeptools/deeptools/parserCommon.py
index 605f22e5c4..5bf6a11399 100644
--- a/pydeeptools/deeptools/parserCommon.py
+++ b/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 '
diff --git a/pydeeptools/deeptools/test/test_bamCoverage_and_bamCompare.py b/pydeeptools/deeptools/test/test_bamCoverage_and_bamCompare.py
index d4ee5a604d..909d513c41 100644
--- a/pydeeptools/deeptools/test/test_bamCoverage_and_bamCompare.py
+++ b/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)
diff --git a/pydeeptools/deeptools/test/test_plotcorrelation.py b/pydeeptools/deeptools/test/test_plotcorrelation.py
index ab30e57018..03dcadd9ee 100644
--- a/pydeeptools/deeptools/test/test_plotcorrelation.py
+++ b/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"
 )
diff --git a/pydeeptools/deeptools/writeBedGraph_bam_and_bw.py b/pydeeptools/deeptools/writeBedGraph_bam_and_bw.py
index 9e47b9fc82..1cb6d4f714 100644
--- a/pydeeptools/deeptools/writeBedGraph_bam_and_bw.py
+++ b/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]))
diff --git a/src/covcalc.rs b/src/covcalc.rs
index a3919a831e..836e424e21 100644
--- a/src/covcalc.rs
+++ b/src/covcalc.rs
@@ -24,7 +24,10 @@ pub fn parse_regions(region: &str, bam_ifile: Vec<&str>) -> (Vec<Region>, HashMa
             .iter()
             .map(|x| {
                 String::from_utf8(x.to_vec()).unwrap_or_else(|e| {
-                    panic!("BAM header for '{}' has a non-UTF-8 chromosome name: {}", bam, e)
+                    panic!(
+                        "BAM header for '{}' has a non-UTF-8 chromosome name: {}",
+                        bam, e
+                    )
                 })
             })
             .collect();
@@ -33,9 +36,9 @@ pub fn parse_regions(region: &str, bam_ifile: Vec<&str>) -> (Vec<Region>, HashMa
             if !found_chroms.contains_key(chrom) {
                 found_chroms.insert(chrom.clone(), 1);
             } else {
-                let count = found_chroms
-                    .get_mut(chrom)
-                    .expect("Chromosome key vanished from found_chroms map between check and update");
+                let count = found_chroms.get_mut(chrom).expect(
+                    "Chromosome key vanished from found_chroms map between check and update",
+                );
                 *count += 1;
             }
         }
@@ -105,9 +108,9 @@ pub fn parse_regions(region: &str, bam_ifile: Vec<&str>) -> (Vec<Region>, HashMa
                 "Supplied chromosome {} is not found.",
                 chromname
             );
-            let chromlen = chromsizes
-                .get(&chromname)
-                .unwrap_or_else(|| panic!("Chromosome '{}' not found in chromsizes map", chromname));
+            let chromlen = chromsizes.get(&chromname).unwrap_or_else(|| {
+                panic!("Chromosome '{}' not found in chromsizes map", chromname)
+            });
             let _reg = Region {
                 chrom: chromname.to_string(),
                 start: Revalue::U(0),
@@ -133,9 +136,9 @@ pub fn parse_regions(region: &str, bam_ifile: Vec<&str>) -> (Vec<Region>, HashMa
                 "Supplied chromosome {} is not found.",
                 chromname
             );
-            let chromlen = chromsizes
-                .get(&chromname)
-                .unwrap_or_else(|| panic!("Chromosome '{}' not found in chromsizes map", chromname));
+            let chromlen = chromsizes.get(&chromname).unwrap_or_else(|| {
+                panic!("Chromosome '{}' not found in chromsizes map", chromname)
+            });
             assert!(
                 end <= *chromlen,
                 "Suplied region end goes beyond chromosome boundary. {} > {}",
@@ -531,7 +534,6 @@ pub fn bam_pileup<'a>(
         // There are two scenarios:
         // bamCoverage mode -> we can collapse bins with same coverage (collapse = true)
         // bamCompare & others -> We cannot collapse the bins, yet. (collapse = false)
-        // Note that collapse can also be passed as a CLI, for those that want that.
         let mut outbuf = Vec::with_capacity(smoothed.len().max(1) * 64);
         let mut push_line = |chrom: &str, s: u32, e: u32, v: f32| {
             use std::io::Write;
@@ -743,12 +745,12 @@ impl Region {
                     anchorstop = *end;
                 }
                 (Revalue::V(start), Revalue::V(end)) => {
-                    anchorstart = *start
-                        .first()
-                        .unwrap_or_else(|| panic!("Region '{}' has an empty exon-start vector", self.name));
-                    anchorstop = *end
-                        .last()
-                        .unwrap_or_else(|| panic!("Region '{}' has an empty exon-end vector", self.name));
+                    anchorstart = *start.first().unwrap_or_else(|| {
+                        panic!("Region '{}' has an empty exon-start vector", self.name)
+                    });
+                    anchorstop = *end.last().unwrap_or_else(|| {
+                        panic!("Region '{}' has an empty exon-end vector", self.name)
+                    });
                 }
                 _ => panic!(
                     "Start and End are not either both u32, or Vecs. This means your regions file is ill-defined. Fix {}.",
@@ -1202,9 +1204,9 @@ impl Region {
                         }
                         if scale_regions.unscaled3prime > 0 {
                             let mut walked_bps: u32 = 0;
-                            let mut lastanchor: u32 = *end
-                                .last()
-                                .unwrap_or_else(|| panic!("Region '{}' has an empty exon-end vector", self.name));
+                            let mut lastanchor: u32 = *end.last().unwrap_or_else(|| {
+                                panic!("Region '{}' has an empty exon-end vector", self.name)
+                            });
                             while walked_bps < scale_regions.unscaled3prime {
                                 let (bin, retanch) = refpoint_exonwalker(
                                     &exons,
@@ -1330,9 +1332,9 @@ impl Region {
 
                         if scale_regions.unscaled5prime > 0 {
                             let mut walked_bps: u32 = 0;
-                            let mut lastanchor: u32 = *end
-                                .last()
-                                .unwrap_or_else(|| panic!("Region '{}' has an empty exon-end vector", self.name));
+                            let mut lastanchor: u32 = *end.last().unwrap_or_else(|| {
+                                panic!("Region '{}' has an empty exon-end vector", self.name)
+                            });
                             while walked_bps < scale_regions.unscaled5prime {
                                 let (bin, retanch) = refpoint_exonwalker(
                                     &exons,
@@ -1757,7 +1759,7 @@ pub struct Gtfparse {
     pub txniddesignator: String,
 }
 
-#[derive(Clone)]
+#[derive(Clone, PartialEq)]
 pub enum Revalue {
     U(u32),
     V(Vec<u32>),
@@ -1832,10 +1834,16 @@ impl Bin {
             Bin::Conbin(start, _) => *start,
             Bin::PaddedConbin(start, _, _) => *start,
             Bin::Catbin(starts) => {
-                starts.first().expect("Bin::Catbin has an empty start vector").0
+                starts
+                    .first()
+                    .expect("Bin::Catbin has an empty start vector")
+                    .0
             }
             Bin::PaddedCatbin(starts, _) => {
-                starts.first().expect("Bin::PaddedCatbin has an empty start vector").0
+                starts
+                    .first()
+                    .expect("Bin::PaddedCatbin has an empty start vector")
+                    .0
             }
         }
     }
@@ -1845,7 +1853,9 @@ impl Bin {
             Bin::PaddedConbin(_, end, _) => *end,
             Bin::Catbin(ends) => ends.last().expect("Bin::Catbin has an empty end vector").1,
             Bin::PaddedCatbin(ends, _) => {
-                ends.last().expect("Bin::PaddedCatbin has an empty end vector").1
+                ends.last()
+                    .expect("Bin::PaddedCatbin has an empty end vector")
+                    .1
             }
         }
     }
diff --git a/src/filehandler.rs b/src/filehandler.rs
index 703b263b2e..45ba0b5b47 100644
--- a/src/filehandler.rs
+++ b/src/filehandler.rs
@@ -373,6 +373,13 @@ pub fn read_bedfile(
                     continue;
                 }
                 end = end.min(chromlen);
+                if start >= end {
+                    println!(
+                        "Warning, region {} has a start position ({}) that is not strictly less than its end position ({}). BED start/end must be 5' -> 3' (start < end); strand belongs in a separate column (use BED6/BED12/GTF to encode strand). Skipping.",
+                        entryname, start, end
+                    );
+                    continue;
+                }
                 if names.contains_key(&entryname) {
                     let count = names
                         .get_mut(&entryname)
@@ -426,6 +433,13 @@ pub fn read_bedfile(
                     continue;
                 }
                 end = end.min(chromlen);
+                if start >= end {
+                    println!(
+                        "Warning, region {} has a start position ({}) that is not strictly less than its end position ({}). BED start/end must be 5' -> 3' (start < end); strand belongs in a separate column (use BED6/BED12/GTF to encode strand). Skipping.",
+                        entryname, start, end
+                    );
+                    continue;
+                }
                 if names.contains_key(&entryname) {
                     let count = names
                         .get_mut(&entryname)
@@ -541,6 +555,13 @@ pub fn read_bedfile(
                             )
                         })
                         .min(chromlen);
+                    if start >= end {
+                        println!(
+                            "Warning, region {} has a start position ({}) that is not strictly less than its end position ({}). BED start/end must be 5' -> 3' (start < end); strand belongs in a separate column. Skipping.",
+                            entryname, start, end
+                        );
+                        continue;
+                    }
                     regions.push(Region {
                         chrom: fields[0].to_string(),  //chrom
                         start: Revalue::U(start),      //start
diff --git a/src/tests/mod.rs b/src/tests/mod.rs
index 8849ebc500..b27ecd53dd 100644
--- a/src/tests/mod.rs
+++ b/src/tests/mod.rs
@@ -1,4 +1,5 @@
 pub mod test_calc;
 pub mod test_covcalc;
+pub mod test_filehandler;
 pub mod test_filtering;
 pub mod test_normalization;
diff --git a/src/tests/test_filehandler.rs b/src/tests/test_filehandler.rs
new file mode 100644
index 0000000000..03db9f38be
--- /dev/null
+++ b/src/tests/test_filehandler.rs
@@ -0,0 +1,60 @@
+use crate::covcalc::Revalue;
+use crate::filehandler::read_bedfile;
+use std::collections::HashMap;
+use std::io::Write;
+use tempfile::NamedTempFile;
+
+fn write_bed(contents: &str) -> NamedTempFile {
+    let mut f = NamedTempFile::new().expect("Failed to create temp BED file");
+    write!(f, "{}", contents).expect("Failed to write temp BED file");
+    f
+}
+
+#[test]
+fn test_read_bedfile_skips_start_ge_end_bed3() {
+    let bed = write_bed("chr1\t100\t200\nchr1\t8000\t3000\nchr1\t300\t400\n");
+    let mut chroms: HashMap<String, u32> = HashMap::new();
+    chroms.insert("chr1".to_string(), 1_000_000);
+
+    let (regions, (_label, entries)) =
+        read_bedfile(&bed.path().to_string_lossy().into_owned(), false, &chroms);
+
+    assert_eq!(entries, 2);
+    assert_eq!(regions.len(), 2);
+    assert_eq!(regions[0].start, Revalue::U(100));
+    assert_eq!(regions[0].end, Revalue::U(200));
+    assert_eq!(regions[0].regionlength, 100);
+    assert_eq!(regions[1].start, Revalue::U(300));
+    assert_eq!(regions[1].end, Revalue::U(400));
+    assert_eq!(regions[1].regionlength, 100);
+}
+
+#[test]
+fn test_read_bedfile_skips_start_ge_end_bed6() {
+    let bed = write_bed(
+        "chr1\t100\t200\tgeneA\t.\t+\nchr1\t13714\t11649\tgeneB\t.\t-\nchr1\t300\t400\tgeneC\t.\t-\n",
+    );
+    let mut chroms: HashMap<String, u32> = HashMap::new();
+    chroms.insert("chr1".to_string(), 1_000_000);
+
+    let (regions, (_label, entries)) =
+        read_bedfile(&bed.path().to_string_lossy().into_owned(), false, &chroms);
+
+    assert_eq!(entries, 2);
+    assert_eq!(regions.len(), 2);
+    assert_eq!(regions[0].name, "geneA");
+    assert_eq!(regions[1].name, "geneC");
+}
+
+#[test]
+fn test_read_bedfile_skips_start_eq_end() {
+    let bed = write_bed("chr1\t100\t200\nchr1\t150\t150\n");
+    let mut chroms: HashMap<String, u32> = HashMap::new();
+    chroms.insert("chr1".to_string(), 1_000_000);
+
+    let (regions, (_label, entries)) =
+        read_bedfile(&bed.path().to_string_lossy().into_owned(), false, &chroms);
+
+    assert_eq!(entries, 1);
+    assert_eq!(regions.len(), 1);
+}