diff --git a/.github/filters.yml b/.github/filters.yml index f703258c3..043116145 100644 --- a/.github/filters.yml +++ b/.github/filters.yml @@ -18,8 +18,8 @@ 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 @@ -27,19 +27,22 @@ deconvolve: &deconvolve - *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,22 +65,28 @@ 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 @@ -83,7 +94,9 @@ naibr: &naibr - '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 @@ -91,11 +104,13 @@ phase: &phase - '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 @@ -103,6 +118,8 @@ assembly: &assembly - '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' diff --git a/harpy/commands/impute.py b/harpy/commands/impute.py index 054222bfc..d6b1c8472 100644 --- a/harpy/commands/impute.py +++ b/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: diff --git a/harpy/commands/phase.py b/harpy/commands/phase.py index 0633b0840..f577d13f0 100644 --- a/harpy/commands/phase.py +++ b/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 diff --git a/harpy/report/notebooks/align_stats.ipynb b/harpy/report/notebooks/align_stats.ipynb index a9fae81a9..9f04250f1 100644 --- a/harpy/report/notebooks/align_stats.ipynb +++ b/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", diff --git a/harpy/validation/vcf.py b/harpy/validation/vcf.py index 5209b928d..fe341144f 100644 --- a/harpy/validation/vcf.py +++ b/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) \ No newline at end of file diff --git a/harpy/validation/xam.py b/harpy/validation/xam.py index 1027404a0..e4db6643c 100644 --- a/harpy/validation/xam.py +++ b/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) + + + \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index e2e6f4492..5d915ea8a 100644 --- a/pyproject.toml +++ b/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"