filter alignments if using --vcf-samples

.github/filters.yml

@@ -18,28 +18,31 @@ validate: &validate
   - 'harpy/snakefiles/validate_bam.smk'
   - 'test/fastq/**'
   - 'test/bam/**'
-  - 'harpy/scripts/check_bam.py'
-  - 'harpy/scripts/check_fastq.py'
+  - 'harpy/utils/check_bam.py'
+  - 'harpy/utils/check_fastq.py'
   - 'harpy/report/notebooks/validate_fastq.ipynb'
   - 'harpy/report/notebooks/validate_bam.ipynb'
 deconvolve: &deconvolve
   - *common
   - *environments
   - 'harpy/commands/deconvolve.py'
   - 'harpy/snakefiles/deconvolve.smk'
+  - 'harpy/validation/fastq.py'
 preproc: &preproc
   - *common
   - *environments
+  - 'haplotag-acbd.py'
   - 'harpy/commands/preprocess.py'
   - 'harpy/snakefiles/preprocess**.smk'
-  - 'haplotag-acbd.py'
+  - 'harpy/validation/fastq.py'
   - 'test/data/preproc/**'
 qc: &qc
   - *common
   - *environments
   - 'harpy/commands/qc.py'
   - 'harpy/snakefiles/qc.smk'
-  - 'harpy/scripts/count_bx.py'
+  - 'harpy/utils/bx_stats_fq.py'
+  - 'harpy/validation/fastq.py'
   - 'test/fastq/**'
 align: &align
   - *common
@@ -50,8 +53,10 @@ align: &align
   - 'harpy/report/notebooks/align_stats.ipynb'
   - 'harpy/report/notebooks/align_lrstats.ipynb'
   - 'harpy/report/notebooks/samtools_stats.ipynb'
-  - 'harpy/scripts/bx_stats.py'
-  - 'harpy/scripts/molecule_coverage.py'
+  - 'harpy/utils/bx_stats_sam.py'
+  - 'harpy/utils/molecule_coverage.py'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/fastq.py'
   - 'test/fastq/**'
 snp: &snp
   - *common
@@ -60,49 +65,61 @@ snp: &snp
   - 'harpy/snakefiles/snp_mpileup.smk'
   - 'harpy/snakefiles/snp_freebayes.smk'
   - 'harpy/report/notebooks/bcftools_stats.ipynb'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/xam.py'
   - 'test/bam/**'
 impute: &impute
   - *common
   - *environments
   - 'harpy/commands/impute.py'
   - 'harpy/snakefiles/impute.smk'
-  - 'test/bam/**'
-  - 'test/vcf/test.bcf'
   - 'harpy/report/notebooks/impute.ipynb'
   - 'harpy/report/notebooks/stitch_collate.ipynb'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/xam.py'
+  - 'test/bam/**'
+  - 'test/vcf/test.bcf'
 leviathan: &leviathan
   - *common
   - *environments
   - 'harpy/commands/sv.py'
   - 'harpy/snakefiles/sv_leviathan.smk'
   - 'harpy/report/notebooks/sv.ipynb'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/xam.py'
   - 'test/bam/**'
 naibr: &naibr
   - *common
   - *environments
   - 'harpy/commands/sv.py'
   - 'harpy/snakefiles/sv_naibr.smk'
   - 'harpy/report/notebooks/sv.ipynb'
-  - 'harpy/scripts/infer-sv.py'
+  - 'harpy/utils/infer_sv.py'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/xam.py'
   - 'test/bam_phased/**'
 phase: &phase
   - *common
   - *environments
   - 'harpy/commands/phase.py'
   - 'harpy/snakefiles/phase_snp.smk'
   - 'harpy/snakefiles/phase_bam.smk'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/xam.py'
   - 'test/bam/**'
   - 'test/vcf/test.bcf'
   - 'test/vcf/test.phased.bcf'
   - 'harpy/report/notebooks/hapcut.ipynb'
-  - 'harpy/scripts/parse_phaseblocks.py'
+  - 'harpy/utils/parse_phaseblocks.py'
 assembly: &assembly
   - *common
   - *environments
   - 'harpy/commands/assembly.py'
   - 'harpy/commands/metassembly.py'
   - 'harpy/snakefiles/assembly.smk'
   - 'harpy/snakefiles/metassembly.smk'
+  - 'harpy/validation/fasta.py'
+  - 'harpy/validation/fastq.py'
 other: &other
   - *common
   - 'harpy/commands/template.py'

harpy/commands/impute.py

@@ -65,6 +65,8 @@ def impute(parameters, vcf, inputs, output, strategy, buffer, grid_size, threads
     vcffile = VCF(vcf, workflow.workflow_directory, quiet)
     vcffile.find_biallelic_contigs()
     vcffile.match_samples(alignments.files, vcf_samples)
+    if vcf_samples:
+        alignments.filter(vcffile.samples)
     region = strategy.lower() != "all" and not strategy.lower().startswith("window:")
     window = None
     if region:

harpy/commands/phase.py

@@ -60,6 +60,8 @@ def bam(vcf, inputs, output, threads, unlinked, vcf_samples, molecule_distance,
     vcffile = VCF(vcf, workflow.workflow_directory, quiet = quiet)
     vcffile.check_phase()
     vcffile.match_samples(alignments.files, vcf_samples)
+    if vcf_samples:
+        alignments.filter(vcffile.samples)
     fasta = FASTA(reference, quiet)
 
     workflow.linkedreads["type"] = alignments.lr_type
@@ -128,6 +130,8 @@ def snp(vcf, inputs, output, threads, unlinked, min_map_quality, min_base_qualit
     vcffile.match_samples(alignments.files, vcf_samples)
     if contigs:
         vcffile.match_contigs(contigs)
+    if vcf_samples:
+        alignments.filter(vcffile.samples)
 
     workflow.linkedreads["type"] = alignments.lr_type
     workflow.notebooks["skip"] = skip_reports

harpy/report/notebooks/align_stats.ipynb

@@ -318,7 +318,7 @@
     "    mask = valids['reads'] >= 2\n",
     "    n_linked_reads = valids.filter(mask)['reads'].sum()\n",
     "    if n_linked_reads == 0:\n",
-    "        print_html(\"No linked molecules (>=2 reads) were detected; skipping linked-only distributions.\")\n",
+    "        print_html(\"No linked molecules (>=2 reads) were detected; linked-read metrics will not be reported.\")\n",
     "        SKIPLR = True\n",
     "    else:\n",
     "        perc_valid = round(tot_valid / n_reads, 4)\n",
@@ -327,13 +327,11 @@
     "\n",
     "    _valid = pl.DataFrame({\"valids\": ['valid', 'invalid'], \"count\": [tot_valid, n_reads - tot_valid]})\n",
     "    _linked_abs = pl.DataFrame({\"linked (absolute)\": ['linked', 'unlinked'], \"count\": [n_linked_reads, n_reads - n_linked_reads]})\n",
-    "    _linked_rel = pl.DataFrame({\"linked (valids)\": ['linked', 'unlinked'], \"count\": [n_linked_reads, tot_valid - n_linked_reads]})\n",
     "\n",
     "    (\n",
     "        alt.hconcat(\n",
     "            piechart(_valid, goodval = 'valid', title = \"Percent Valid\"),\n",
-    "            piechart(_linked_abs, goodval = \"linked\", title =  \"Linked (all)\"),\n",
-    "            piechart(_linked_rel, goodval = \"linked\", title =  \"Linked (valid)\"),\n",
+    "            piechart(_linked_abs, goodval = \"linked\", title =  \"Percent Linked\"),\n",
     "            convolutionbar(pl.DataFrame(_d), title = \"Deconvolutions\"),\n",
     "            center = True,\n",
     "            bounds = 'flush'\n",

harpy/validation/vcf.py

@@ -9,7 +9,6 @@
 
 from harpy.common.printing import HarpyPrint
 
-
 class VCF():
     '''
     A class to contain and validate a VCF input file.
@@ -23,6 +22,11 @@ def __init__(self, filename:str, workdir:str, quiet:int = 0):
         self.contigs: dict[str,int] = {}
         self.print = HarpyPrint(quiet)
 
+        if self.file.lower().endswith("bcf") and not os.path.exists(f"{self.file}.csi"):
+            pysam.bcftools.index(self.file)
+        if self.file.lower().endswith("vcf.gz") and not os.path.exists(f"{self.file}.tbi"):
+            pysam.bcftools.index("--tbi", self.file)
+
     def get_contigs(self):
         """reads the header of a vcf/bcf file and populate `self.contigs` with the contigs (keys) and their lengths (values)"""
         with pysam.VariantFile(self.file) as _vcf:
@@ -39,25 +43,19 @@ def find_biallelic_contigs(self):
         self.biallelic_file = Path(os.path.join(self.workdir, os.path.basename(self.file) + ".biallelic")).resolve().as_posix()
         if not self.contigs:
             self.get_contigs()
-        to_keep = []
-
-        if self.file.lower().endswith("bcf") and not os.path.exists(f"{self.file}.csi"):
-            pysam.bcftools.index(self.file)
-        if self.file.lower().endswith("vcf.gz") and not os.path.exists(f"{self.file}.tbi"):
-            pysam.bcftools.index("--tbi", self.file)
 
         bcftools = which("bcftools") or "bcftools"
 
         with open(self.biallelic_file, "w", encoding="utf-8") as f:
             for contig in list(self.contigs.keys()):
+                snpcount = 0
+                keep = False
                 # Use bcftools to count the number of biallelic SNPs in the contig
                 viewcmd = subprocess.Popen(
                     [bcftools, 'view', '-H', '-r', str(contig), '-v', 'snps', '-m2', '-M2', '-c', '2', self.file],
                     stdout=subprocess.PIPE,
                     text = True
                 )
-                snpcount = 0
-                keep = False
                 while True:
                     # Read the next line of output
                     line = viewcmd.stdout.readline()
@@ -181,3 +179,8 @@ def validate_region(self, region: str) -> tuple[str,int,int]:
         else:
             self.print.validation(True)
         return contig, startpos, endpos
+
+    def samplelist(self) -> list[str]:
+        '''return the list of sample names from the vcf file'''
+        with pysam.VariantFile(self.file) as _vcf:
+            return list(_vcf.header.samples)

harpy/validation/xam.py

@@ -153,3 +153,16 @@ def which_linkedread(self, file_path: str) -> str:
                 if HAPLOTAGGING_RX.search(bx):
                     return "haplotagging"
         return "none"
+
+    def filter(self, keeplist: list[str]):
+        '''Filter self.files by their basenames to keep only those samples in `keeplist`'''
+        tmp = []
+        re_ext = re.compile(r"\.(bam|sam)$", re.IGNORECASE)
+        for i in self.files:
+            if os.path.basename(re_ext.sub("", str(i))) in keeplist:
+                tmp.append(i)
+        self.files = tmp
+        self.count = len(tmp)
+
+
+

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "harpy"
-version = "4.1.4"
+version = "4.1.5"
 authors = [{name = "Pavel Dimens"}]
 description = "Automated workflows to preprocess sequences, trim reads, map sequences, call variants, impute genotypes, and phase haplotypes of Haplotagging data. Batteries included."
 requires-python = ">=3.11"