filter alignments if using --vcf-samples

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/validation/vcf.py

@@ -39,7 +39,6 @@ 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)
@@ -181,3 +180,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,15 @@ 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
+
+
+