gz bed/gtf and test

Author: WardDeb

pydeeptools/deeptools/test/test_computeMatrix_referencepoint.py

@@ -13,6 +13,8 @@
 REGIONS_IN1 = ROOT + "input_computeMatrix_regions1.bed"
 REGIONS_IN2 = ROOT + "input_computeMatrix_regions2.bed"
 REGIONS_GTF = ROOT + "input_computeMatrix_regions3.gtf"
+REGIONS_IN1_GZ = ROOT + "input_computeMatrix_regions1.bed.gz"
+REGIONS_GTF_GZ = ROOT + "input_computeMatrix_regions3.gtf.gz"
 REGIONS_BED12 = ROOT + "input_computeMatrix_regions4bed12.bed"
 BL = ROOT + "input_computeMatrix_blacklist.bed"
 BIGWIG_IN1 = ROOT + "input_computeMatrix_bw1.bw"
@@ -246,6 +248,28 @@ def test_compute_matrix_refpoint():
     mat_match, mat_diffs = _compare_tab_files(outfile_mat, exp_mat)
     assert mat_match, f"matrix mismatch: {mat_diffs}"
 
+def test_compute_matrix_refpoint_bed_gz():
+    exp_npz = ROOT + "mat.gz"
+    exp_mat = ROOT + "mat.tab"
+
+    _, outfile_npz = tempfile.mkstemp(suffix='.gz')
+    _, outfile_mat = tempfile.mkstemp(suffix='.tab')
+    args = f"""
+    reference-point
+    --regionsFileName {REGIONS_IN1_GZ}
+    --scoreFileName {BIGWIG_IN1}
+    -b 20 -a 20
+    -o {outfile_npz} --outFileNameMatrix {outfile_mat}
+    """.split()
+    cm.main(args)
+
+    header_differences, data_differences, rowdiffdic = _compare_mat_gz(outfile_npz, exp_npz)
+    assert not header_differences, f"header mismatch: {header_differences}"
+    assert not data_differences, f"data mismatch: {data_differences}\nrowdiffdict: {rowdiffdic}"
+
+    mat_match, mat_diffs = _compare_tab_files(outfile_mat, exp_mat)
+    assert mat_match, f"matrix mismatch: {mat_diffs}"
+
 def test_compute_matrix_refpoint2():
     exp_npz = ROOT + "mat2.gz"
     exp_mat = ROOT + "mat2.tab"
@@ -958,6 +982,28 @@ def test_compute_matrix_refpoint_gtf():
     mat_match, mat_diffs = _compare_tab_files(outfile_mat, exp_mat)
     assert mat_match, f"matrix mismatch: {mat_diffs}"
 
+def test_compute_matrix_refpoint_gtf_gz():
+    exp_npz = ROOT + "mat_gtf.gz"
+    exp_mat = ROOT + "mat_gtf.tab"
+
+    _, outfile_npz = tempfile.mkstemp(suffix='.gz')
+    _, outfile_mat = tempfile.mkstemp(suffix='.tab')
+    args = f"""
+    reference-point
+    --regionsFileName {REGIONS_GTF_GZ}
+    --scoreFileName {BIGWIG_IN1} {BIGWIG_IN2} {BIGWIG_IN3}
+    -b 30 -a 30 --binSize 30
+    -o {outfile_npz} --outFileNameMatrix {outfile_mat}
+    """.split()
+    cm.main(args)
+
+    header_differences, data_differences, rowdiffdic = _compare_mat_gz(outfile_npz, exp_npz)
+    assert not header_differences, f"header mismatch: {header_differences}"
+    assert not data_differences, f"data mismatch: {data_differences}\nrowdiffdict: {rowdiffdic}"
+
+    mat_match, mat_diffs = _compare_tab_files(outfile_mat, exp_mat)
+    assert mat_match, f"matrix mismatch: {mat_diffs}"
+
 def test_compute_matrix_refpoint_bed12mg():
     exp_npz = ROOT + "mat_bed12mg.gz"
     exp_mat = ROOT + "mat_bed12mg.tab"

pydeeptools/deeptools/test/test_data/test_computematrix/input_computeMatrix_regions1.bed.gz

@@ -18,6 +18,8 @@
 BAMB = ROOT + "testB"
 BLACKLIST = ROOT + "blacklist.bed"
 GTF2 = ROOT + "test_customids.gtf"
+BED_GZ = ROOT + "test.bed.gz"
+GTF_GZ = ROOT + "test.gtf.gz"
 
 def compare_tsv(exp_tsv, obs_tsv, delta=1.0):
     with open(exp_tsv) as e, open(obs_tsv) as o:
@@ -63,6 +65,32 @@ def test_multiBamSummary_bedmode_bed():
             )
         )
 
+def test_multiBamSummary_bedmode_gtf_gz():
+    _, outfile = tempfile.mkstemp(suffix=".npz")
+    for fname in ['.bam', '.cram']:
+        fname = BAM + fname
+        args = f"BED-file --BED {GTF_GZ} -b {fname} {fname} -o {outfile}".split()
+        mbs.main(args)
+        resp = np.load(outfile)
+        matrix = resp['matrix']
+
+        nt.assert_allclose(matrix, np.array([[144.0, 144.0],[143.0, 143.0]]))
+
+def test_multiBamSummary_bedmode_bed_gz():
+    _, outfile = tempfile.mkstemp(suffix=".npz")
+    for fname in ['.bam', '.cram']:
+        fname = BAM + fname
+        args = f"BED-file --BED {BED_GZ} -b {fname} {fname} -o {outfile}".split()
+        mbs.main(args)
+        resp = np.load(outfile)
+        matrix = resp['matrix']
+
+        nt.assert_allclose(matrix,
+            np.array(
+                [[1.0, 1.0], [144.0, 144.0], [144.0, 144.0], [6.0, 6.0], [143.0, 143.0], [22.0, 22.0], [25.0, 25.0], [1.0, 1.0], [0.0, 0.0]]
+            )
+        )
+
 def test_multiBamSummary_bedmode_multibed():
     _, outfile = tempfile.mkstemp(suffix=".npz")
     for fname in ['.bam', '.cram']:

pydeeptools/deeptools/test/test_data/test_computematrix/input_computeMatrix_regions3.gtf.gz

@@ -3,6 +3,7 @@ use crate::covcalc::{Bin, Gtfparse, Region, Revalue, Scalingregions};
 use bigtools::beddata::BedParserStreamingIterator;
 use bigtools::{BigWigRead, BigWigWrite, Value};
 use flate2::Compression;
+use flate2::read::MultiGzDecoder;
 use flate2::write::GzEncoder;
 use itertools::Itertools;
 use rust_htslib::bam::{IndexedReader, Read, Reader};
@@ -55,10 +56,22 @@ where
     }
 }
 
+fn open_bed_or_gtf_reader(fp: &str) -> BufReader<Box<dyn std::io::Read>> {
+    let mut file = File::open(fp).expect(format!("Failed to open file: {}", fp).as_str());
+    let mut magic = [0u8; 2];
+    let n = std::io::Read::read(&mut file, &mut magic).unwrap_or(0);
+    std::io::Seek::seek(&mut file, std::io::SeekFrom::Start(0))
+        .expect(format!("Failed to seek in file: {}", fp).as_str());
+    if n == 2 && magic[0] == 0x1f && magic[1] == 0x8b {
+        BufReader::new(Box::new(MultiGzDecoder::new(file)) as Box<dyn std::io::Read>)
+    } else {
+        BufReader::new(Box::new(file) as Box<dyn std::io::Read>)
+    }
+}
+
 pub fn is_bed_or_gtf(fp: &str) -> String {
     // Check if file is a bed or gtf file.
-    let file = File::open(fp).expect(format!("Failed to open file: {}", fp).as_str());
-    let reader = BufReader::new(file);
+    let reader = open_bed_or_gtf_reader(fp);
     // Get the first line that doesn't start with #
     for line in reader.lines() {
         let line = line.unwrap();
@@ -85,7 +98,7 @@ pub fn read_gtffile(
     let mut entries: u32 = 0;
     let mut txnids: Vec<String> = Vec::new();
 
-    let gtffile = BufReader::new(File::open(gtf_file).unwrap());
+    let gtffile = open_bed_or_gtf_reader(gtf_file);
 
     if gtfparse.metagene {
         // metagene implementation - more work here.
@@ -157,9 +170,6 @@ pub fn read_gtffile(
             let has_exons = txn_hash.contains_key(&txnid);
             let has_transcript = txn_transcript.contains_key(&txnid);
 
-            // Only keep entries that have a 'transcript' feature line.
-            // This mirrors Python behavior: the transcript must exist in the
-            // interval tree for exons to be returned by findOverlaps().
             if has_exons && !has_transcript {
                 continue;
             }
@@ -266,7 +276,7 @@ pub fn read_bedfile(
     let mut nonbed12: bool = false;
     let mut entries: u32 = 0;
 
-    let bedfile = BufReader::new(File::open(bed_file).unwrap());
+    let bedfile = open_bed_or_gtf_reader(bed_file);
 
     for line in bedfile.lines() {
         let line = line.unwrap();
@@ -1011,7 +1021,6 @@ pub fn write_sorted_regions_bed(
         group_boundaries.push(cumsum);
     }
     for (idx, region) in regions.iter().enumerate() {
-        // Find label_idx: last boundary <= idx, matching Python: np.flatnonzero(boundaries <= idx)[-1]
         let label_idx = group_boundaries
             .iter()
             .take_while(|&&b| b <= idx as u32)