ci: improve test coverage

Author: js2264

.gitignore

@@ -64,7 +64,6 @@ tmp/
 tests_data/tmp
 tests_data/out_test
 tests_data/outdir/plot/
-tests_data/outdir_pipeline/
 data/
 
 # Translations

src/metator/commands.py

@@ -1011,7 +1011,9 @@ class Qc(AbstractCommand):
     usage:
         qc --assembly=FILE --enzyme=STR [--bin-summary=FILE]
         [--contig-data=FILE] [--metator-dir=DIR] [--no-clean-up] [--outdir=DIR]
-        [--prefix=STR] [--plot] [--threshold=STR] [--tmpdir=DIR] <pairsfile>...
+        [--prefix=STR] [--plot] [--threshold=STR] [--tmpdir=DIR]
+        [--completeness=FLOAT] [--redundancy=FLOAT]
+        <pairsfile>...
 
     arguments:
         pairsfile           File(s) containing pairs information.
@@ -1022,6 +1024,8 @@ class Qc(AbstractCommand):
                                 contigs of one bin.
         -b, --bin-summary=FILE  Path to the bin_summary.txt file from MetaTOR
                                 output.
+        -C, --completeness=FLOAT  Minimum completeness of the bins to consider
+                                them as high quality. [Default: 0.7]
         -c, --contig-data=FILE  Path to the contig_data_final.txt file from
                                 MetaTOR output.
         -e, --enzyme=STR        The list of restriction enzyme used to digest
@@ -1034,6 +1038,8 @@ class Qc(AbstractCommand):
         -o, --outdir=DIR        Directory to save output plots and log.
         -p, --prefix=STR        Name of the sample to add on plot and files.
         -P, --plot              If enable display some plots.
+        -R, --redundancy=FLOAT  Minimum redundancy of the bins to consider
+                                them as high quality. [Default: 1.15]
         -t, --threshold=STR     Hicstuff religation and loop thresholds.
                                 Two integers seperated by a coma.
         -T, --tmpdir=DIR        Temporary directory to save pairs files
@@ -1095,6 +1101,14 @@ def execute(self):
             self.args["--enzyme"] = self.args["--enzyme"].split(",")
         if self.args["--threshold"]:
             self.args["--threshold"] = self.args["--threshold"].split(",")
+        if len(self.args["--completeness"]) == 0:
+            self.args["--completeness"] = 0.7
+        else:
+            self.args["--completeness"] = float(self.args["--completeness"])
+        if len(self.args["--redundancy"]) == 0:
+            self.args["--redundancy"] = 1.15
+        else:
+            self.args["--redundancy"] = float(self.args["--redundancy"])
 
         # Launch quality check
         mtq.quality_check(
@@ -1108,6 +1122,8 @@ def execute(self):
             self.args["--plot"],
             enzyme=self.args["--enzyme"],
             threshold=self.args["--threshold"],
+            completeness_threshold=self.args["--completeness"],
+            redundancy_threshold=self.args["--redundancy"],
         )
 
         # Delete the temporary folder.

src/metator/quality_check.py

@@ -3,8 +3,8 @@
 
 """Generates quality metrics from the output of MetaTOR.
 
-General utilities functions to extract high  quality contigs which will be used 
-to estimate the quality of the metaHiC librairies nad generates some plots to 
+General utilities functions to extract high  quality contigs which will be used
+to estimate the quality of the metaHiC librairies nad generates some plots to
 illustrate it.
 
 Core functions to assess the quality:
@@ -26,7 +26,7 @@
 from os.path import join
 
 
-def extract_hq_contigs(bin_summary, contigs_data):
+def extract_hq_contigs(bin_summary, contigs_data, completeness_threshold=0.7, redundancy_threshold=1.15):
     """Function to extract the high quality contigs from the metator output.
     These contigs will be the one used to assess the quality of the dataset.
 
@@ -47,8 +47,8 @@ def extract_hq_contigs(bin_summary, contigs_data):
     # Extract high quality MAGs.
     hq_mags = bin_summary.index[
         np.logical_and(
-            (bin_summary["Weighted completeness"] > 0.7),
-            (bin_summary["Weighted redundancy"] < 1.15),
+            (bin_summary["Weighted completeness"] > completeness_threshold),
+            (bin_summary["Weighted redundancy"] < redundancy_threshold),
         )
     ]
     n_mags = len(hq_mags)
@@ -83,13 +83,12 @@ def extract_pairs(pairs_files, out_file, contigs, contigs_data):
     with open(out_file, "w") as output_pairs:
         # Write the header of the output pairs
         output_pairs.write("## pairs format v1.0\n")
-        output_pairs.write(
-            "#columns: readID chr1 pos1 chr2 pos2 strand1 strand2\n"
-        )
+        output_pairs.write("#columns: readID chr1 pos1 chr2 pos2 strand1 strand2\n")
         for contig in contigs:
             output_pairs.write(
                 "#chromsize: {0} {1}\n".format(
-                    contig, contigs_data.loc[contig, "Size"],
+                    contig,
+                    contigs_data.loc[contig, "Size"],
                 )
             )
         for pairs_file in pairs_files:
@@ -166,9 +165,7 @@ def hic_quality(
     restrict_table = {}
     for record in SeqIO.parse(mio.read_compressed(fasta), "fasta"):
         # Get chromosome restriction table
-        restrict_table[record.id] = hcd.get_restriction_table(
-            record.seq, enzyme, circular=False
-        )
+        restrict_table[record.id] = hcd.get_restriction_table(record.seq, enzyme, circular=False)
 
     # Add fragment index to pairs (readID, chr1, pos1, chr2,
     # pos2, strand1, strand2, frag1, frag2)
@@ -187,11 +184,7 @@ def hic_quality(
             fig_path=plot_event,
             prefix=prefix,
         )
-        logger.info(
-            "Filtering with thresholds: uncuts={0} loops={1}".format(
-                uncut_thr, loop_thr
-            )
-        )
+        logger.info("Filtering with thresholds: uncuts={0} loops={1}".format(uncut_thr, loop_thr))
     # Filter reads and save metrics on informative reads
     n_religated = 0
     n_loops = 0
@@ -230,11 +223,7 @@ def hic_quality(
                 else:
                     n_inter_mags += 1
     if n_intra_mags + n_inter_mags > 0:
-        rat_info = (
-            100
-            * (n_informative_intra + n_informative_inter)
-            / (n_intra_mags + n_inter_mags)
-        )
+        rat_info = 100 * (n_informative_intra + n_informative_inter) / (n_intra_mags + n_inter_mags)
         noise_ratio = 100 * n_inter_mags / (n_inter_mags + n_intra_mags)
     else:
         logger.warning("No pairs have benn extracted. All scores set to 0.")
@@ -251,9 +240,7 @@ def hic_quality(
     logger.info(f"Loop ratio: {100 * n_loops / n_intra_mags:.2f}%.")
     logger.info(f"Weirds ratio: {100 * n_weirds / n_intra_mags:.2f}%.")
     logger.info(f"Informative contacts estimation: {rat_info:.2f}%.")
-    logger.info(
-        f"Ratio inter/intra contigs: {n_informative_inter / (n_informative_intra + n_informative_inter):.2f}%."
-    )
+    logger.info(f"Ratio inter/intra contigs: {n_informative_inter / (n_informative_intra + n_informative_inter):.2f}%.")
     logger.info(f"Noise contact ratio: {noise_ratio:.2f}%")
     logger.info(f"Noise score: {noise_score:.2E}")
 
@@ -293,6 +280,8 @@ def quality_check(
     plot,
     enzyme,
     threshold,
+    completeness_threshold=0.7,
+    redundancy_threshold=1.15,
 ):
     """Main function to compute the quality of the metaHiC library and to
     display some metrics about it.
@@ -327,11 +316,11 @@ def quality_check(
     pairs_idx = join(tmp_dir, f"{prefix}_idx.pairs")
 
     # Extract high quality contigs.
-    hq_contigs, n_mags = extract_hq_contigs(bin_summary, contigs_data)
-    # Extract pairs
-    n_pairs = extract_pairs(
-        pairs_files, pairs, list(hq_contigs.keys()), contigs_data
+    hq_contigs, n_mags = extract_hq_contigs(
+        bin_summary, contigs_data, completeness_threshold=completeness_threshold, redundancy_threshold=redundancy_threshold
     )
+    # Extract pairs
+    n_pairs = extract_pairs(pairs_files, pairs, list(hq_contigs.keys()), contigs_data)
 
     # Estimate HiC quality.
     if n_pairs > 0:

tests/test_commands.py

@@ -33,6 +33,8 @@ def tmp_dir(tmp_path_factory):
     "SUBSET_FASTA": "tests_data/subset.assembly.fa",
     "SUBSET_PAIRS": "tests_data/subset.pairs.gz",
     "SUBSET_MGES": "tests_data/subset.mges.txt",
+    "OUT_PIPELINE": "tests_data/outdir_pipeline",
+    "HOSTMAG_FASTA": "tests_data/outdir_pipeline/metator_00005_00000.fa",
 }
 NORMALIZE = (
     "norm",
@@ -169,6 +171,50 @@ def test_pipeline(tmp_path):
     proc.execute()
 
 
+def test_qc(tmp_path):
+    """Test the metator qc command with subset data."""
+    args = (
+        "--assembly {SUBSET_FASTA} --enzyme HindIII,DpnII --bin-summary {OUT_PIPELINE}/bin_summary.txt --contig-data {OUT_PIPELINE}/contig_data_final.txt  -o {OUT_TEST} -P -R 1.5 {SUBSET_PAIRS}"
+    ).format(
+        OUT_TEST=Path(tmp_path, "out_test"),
+        **global_args,
+    )
+    print("CLI command:\n", "metator qc", *args.split(" "))
+
+    proc = mtc.Qc(args.split(" "), {})
+    proc.execute()
+
+
+def test_contactmap(tmp_path):
+    """Test the metator contactmap command with subset data."""
+    os.makedirs(Path(tmp_path, "out_test"), exist_ok=True)
+    args = (
+        "--assembly {SUBSET_FASTA} --enzyme HindIII,DpnII --contig-data {OUT_PIPELINE}/contig_data_final.txt --name metator_00005_00000 --filter --mat-fmt cool -o {OUT_TEST} --pcr-dup {SUBSET_PAIRS}"
+    ).format(
+        OUT_TEST=Path(tmp_path, "out_test"),
+        **global_args,
+    )
+    print("CLI command:\n", "metator contactmap", *args.split(" "))
+
+    proc = mtc.Contactmap(args.split(" "), {})
+    proc.execute()
+
+
+def test_scaffold(tmp_path):
+    """Test the metator scaffold command with subset data."""
+    os.makedirs(Path(tmp_path, "out_test"), exist_ok=True)
+    args = (
+        "--bin-name metator_00005_00000 --input-fasta {HOSTMAG_FASTA} --out-fasta {OUT_TEST}/metator_00005_00000_scaffolded.fa --out-frags {OUT_TEST}/metator_00005_00000_scaffolded.frags {SUBSET_PAIRS}"
+    ).format(
+        OUT_TEST=Path(tmp_path, "out_test"),
+        **global_args,
+    )
+    print("CLI command:\n", "metator scaffold", *args.split(" "))
+
+    proc = mtc.Scaffold(args.split(" "), {})
+    proc.execute()
+
+
 def test_pairs(tmp_path):
     OUT_TEST = Path(tmp_path, "out_test")
     os.makedirs(OUT_TEST, exist_ok=True)

tests_data/outdir/bin_summary.txt

@@ -1,4 +1,4 @@
-	Length	GC-content	Present	Markers	Completeness	Redundancy	Weighted	completeness	Weighted	redundancy	Contigs	N50	L50	N90	L90	CDs	HiC_abundance
+	Length	GC-content	Present Markers	Completeness	Redundancy	Weighted completeness	Weighted redundancy	Contigs	N50	L50	N90	L90	CDs	HiC_abundance
 MetaTOR_66_0	1513026.00	33.62	85.00	0.81	1.18	0.72	1.24	1656.00	1135.00	388.00	434.00	1234.00	2653.00	0.0597
 MetaTOR_37_0	1681556.00	38.44	77.00	0.73	1.13	0.67	1.18	1364.00	1624.00	334.00	619.00	990.00	2600.00	0.0075
 MetaTOR_14_0	1543324.00	43.89	41.00	0.39	1.20	0.55	1.12	2312.00	742.00	676.00	367.00	1846.00	3223.00	0.0357

tests_data/outdir_pipeline/bin_summary.txt

@@ -0,0 +1,41 @@
+	Length	GC-content	Present Markers	Completeness	Redundancy	Weighted completeness	Weighted redundancy	Contigs	N50	L50	N90	L90	CDs	recursive	step	parent	HiC_abundance
+metator_00041_00000	51418.00	40.32	2.00	0.02	1.00	0.04	1.00	32.00	2090.00	8.00	733.00	24.00	67.00	False	0		0.0190
+metator_00039_00000	35750.00	39.16	0.00	0.00	0.00	0.00	0.00	12.00	4219.00	2.00	1057.00	8.00	40.00	False	0		0.0075
+metator_00003_00000	57791.00	38.04	0.00	0.00	0.00	0.00	0.00	27.00	4565.00	5.00	1062.00	19.00	69.00	False	0		0.0158
+metator_00101_00000	118235.00	37.56	5.00	0.05	1.00	0.05	1.00	5.00	43774.00	2.00	11695.00	3.00	131.00	False	0		0.0042
+metator_00011_00000	42366.00	40.36	2.00	0.02	1.00	0.02	1.00	20.00	3018.00	5.00	938.00	15.00	46.00	False	0		0.0131
+metator_00053_00000	47335.00	39.68	1.00	0.01	1.00	0.02	1.00	20.00	5205.00	2.00	844.00	14.00	60.00	False	0		0.0147
+metator_00002_00000	3361245.00	38.39	104.00	0.99	1.84	0.99	1.51	687.00	13185.00	72.00	1950.00	325.00	4026.00	False	0		21.6089
+metator_00045_00000	62395.00	41.86	0.00	0.00	0.00	0.00	0.00	38.00	2657.00	8.00	723.00	28.00	80.00	False	0		0.0279
+metator_00063_00000	54066.00	41.69	5.00	0.05	1.40	0.03	1.29	20.00	5567.00	2.00	880.00	13.00	76.00	False	0		0.0102
+metator_00069_00000	38265.00	41.83	8.00	0.06	1.12	0.05	1.00	26.00	2232.00	7.00	588.00	20.00	56.00	False	0		0.0146
+metator_00094_00000	37829.00	41.54	0.00	0.00	0.00	0.00	0.00	21.00	2305.00	7.00	992.00	17.00	50.00	False	0		0.0175
+metator_00047_00000	72931.00	39.93	3.00	0.03	1.00	0.01	1.00	29.00	3135.00	5.00	1110.00	19.00	109.00	False	0		0.0354
+metator_00082_00000	32506.00	43.16	0.00	0.00	0.00	0.00	0.00	13.00	3602.00	3.00	876.00	9.00	45.00	False	0		0.0142
+metator_00074_00000	30696.00	42.99	4.00	0.04	1.00	0.04	1.00	23.00	1550.00	8.00	797.00	19.00	52.00	False	0		0.0113
+metator_00022_00000	59382.00	42.07	1.00	0.01	1.00	0.01	1.00	32.00	4774.00	5.00	677.00	23.00	71.00	False	0		0.0237
+metator_00034_00000	89318.00	41.91	6.00	0.06	1.00	0.05	1.00	41.00	4066.00	6.00	773.00	28.00	119.00	False	0		0.0328
+metator_00026_00000	43467.00	42.10	0.00	0.00	0.00	0.00	0.00	29.00	1721.00	10.00	915.00	23.00	62.00	False	0		0.0187
+metator_00040_00000	51486.00	42.38	2.00	0.02	1.00	0.04	1.00	19.00	3870.00	4.00	1075.00	13.00	59.00	False	0		0.0220
+metator_00060_00000	52506.00	41.35	1.00	0.01	1.00	0.01	1.00	24.00	2891.00	6.00	1090.00	18.00	59.00	False	0		0.0131
+metator_00033_00000	31153.00	39.08	2.00	0.02	1.00	0.01	1.00	17.00	2523.00	5.00	829.00	13.00	46.00	False	0		0.0092
+metator_00048_00000	84798.00	41.80	8.00	0.08	1.00	0.10	1.00	56.00	2025.00	16.00	765.00	44.00	120.00	False	0		0.0338
+metator_00078_00000	49578.00	40.35	2.00	0.02	1.50	0.00	1.28	19.00	4314.00	3.00	898.00	14.00	56.00	False	0		0.0145
+metator_00032_00000	52425.00	37.70	2.00	0.02	2.00	0.04	2.00	25.00	2781.00	5.00	872.00	19.00	79.00	False	0		0.0212
+metator_00004_00000	35226.00	40.71	0.00	0.00	0.00	0.00	0.00	21.00	1921.00	5.00	917.00	17.00	54.00	False	0		0.0163
+metator_00005_00000	6815157.00	49.40	88.00	0.84	1.39	0.92	1.36	269.00	57609.00	30.00	13183.00	116.00	6263.00	False	0		27.0418
+metator_00062_00000	37422.00	39.93	3.00	0.02	1.00	0.01	1.00	16.00	4001.00	3.00	844.00	11.00	47.00	False	0		0.0153
+metator_00020_00000	40158.00	42.38	3.00	0.03	1.00	0.01	1.00	24.00	2145.00	6.00	910.00	19.00	56.00	False	0		0.0136
+metator_00100_00000	36558.00	42.09	2.00	0.02	1.50	0.02	1.58	13.00	5301.00	2.00	1039.00	8.00	46.00	False	0		0.0100
+metator_00057_00000	31875.00	43.57	3.00	0.03	1.00	0.03	1.00	25.00	1269.00	8.00	701.00	20.00	49.00	False	0		0.0127
+metator_00006_00000	58284.00	42.92	2.00	0.02	1.00	0.03	1.00	24.00	5559.00	3.00	1015.00	16.00	64.00	False	0		0.0161
+metator_00037_00000	33618.00	41.10	1.00	0.01	1.00	0.02	1.00	17.00	2289.00	3.00	750.00	13.00	47.00	False	0		0.0094
+metator_00066_00000	35002.00	37.61	2.00	0.02	1.00	0.02	1.00	14.00	6705.00	3.00	811.00	9.00	48.00	False	0		0.0112
+metator_00016_00000	80386.00	40.03	1.00	0.01	1.00	0.01	1.00	27.00	7157.00	4.00	1004.00	16.00	102.00	False	0		0.0401
+metator_00007_00000	45121.00	41.47	2.00	0.02	1.00	0.05	1.00	26.00	2420.00	6.00	870.00	20.00	68.00	False	0		0.0201
+metator_00064_00000	33279.00	42.00	0.00	0.00	0.00	0.00	0.00	13.00	4484.00	4.00	1175.00	9.00	45.00	False	0		0.0106
+metator_00049_00000	54634.00	41.81	1.00	0.01	1.00	0.02	1.00	34.00	2165.00	9.00	745.00	26.00	78.00	False	0		0.0249
+metator_00001_00000	49671.00	40.74	2.00	0.02	1.00	0.00	1.00	26.00	2596.00	7.00	939.00	20.00	59.00	False	0		0.0194
+metator_00009_00000	31680.00	42.26	2.00	0.02	1.00	0.01	1.00	14.00	4609.00	3.00	1018.00	9.00	41.00	False	0		0.0091
+metator_00073_00000	34495.00	37.80	2.00	0.02	1.00	0.04	1.00	21.00	2137.00	4.00	699.00	16.00	56.00	False	0		0.0156
+metator_00008_00000	63542.00	41.12	4.00	0.04	1.00	0.04	1.00	42.00	1691.00	14.00	930.00	33.00	97.00	False	0		0.0293