feat: added wrapper for Paraphase

bio/paraphase/environment.yaml

@@ -0,0 +1,10 @@
+channels:
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - minimap2=2.28
+  - paraphase=3.1.1
+  - htslib =1.20
+  - samtools =1.20
+  - bcftools =1.20

bio/paraphase/meta.yaml

@@ -0,0 +1,16 @@
+name: "paraphase"
+name: paraphase
+url: https://github.com/PacificBiosciences/paraphase
+description: paraphase resolves SNVs in gene families for PacBio longreads from whole genomes
+authors: Magdalena
+input:
+  bam: Path to bam file with aligned genomic reads
+  fasta: genome fasta file for target genome
+  faidx: fasta .fai index for the genome fasta file
+output:
+  vcf: vcf file with resulting phased variants
+  bam: bam file with realigned reads
+  bai: bai index file for the realigned reads
+  merged_vcf: VCF file containing all newly discovered variants, phased by gene copy and haplotype
+params:
+  extra: additional parameters that will be forwarded to paraphase

bio/paraphase/test/Snakefile

@@ -0,0 +1,24 @@
+# Set the custom wrapper prefix for the entire workflow
+
+
+rule paraphase:
+    input:
+        bam="test_data/PDPK1.bam",
+        fasta="test_data/chr16s.fa",
+        faidx="test_data/homo_sapiens.fa.fai",
+    output:
+        bam="paraphase/PDPK1_realigned.paraphase.bam",
+        bai="paraphase/PDPK1_realigned.paraphase.bam.bai",
+        merged_vcf="paraphase/PDK1_paraphase.vcf.gz",
+        vcf_header="paraphase/vcf_chromosome_header.vcf",
+    params:
+        extra="-g PDPK1 --genome 38",
+    log:
+        "paraphase/PDPK1_realigned.vcf.log",
+    benchmark:
+        "paraphase/PDPK1_paraphase.vcf.benchmark.tsv"
+    threads: 1
+    resources:
+        mem_mb=1024,
+    wrapper:
+        "main/bio/paraphase"

bio/paraphase/test/test_data/homo_sapiens.fa.fai

@@ -0,0 +1,25 @@
+chr1	2012	6	2012	2013
+chr2	2012	2025	2012	2013
+chr3	2012	4044	2012	2013
+chr4	2012	6063	2012	2013
+chr5	71099999	8082	71099999	71100000
+chr6	2012	71108088	2012	2013
+chr7	2012	71110107	2012	2013
+chr8	2012	71112126	2012	2013
+chr9	2012	71114145	2012	2013
+chr10	2012	71116165	2012	2013
+chr11	2012	71118185	2012	2013
+chr12	2012	71120205	2012	2013
+chr13	2012	71122225	2012	2013
+chr14	2012	71124245	2012	2013
+chr15	2012	71126265	2012	2013
+chr16	2012	71128285	2012	2013
+chr17	2012	71130305	2012	2013
+chr18	2012	71132325	2012	2013
+chr19	2012	71134345	2012	2013
+chr20	2012	71136365	2012	2013
+chr21	2012	71138385	2012	2013
+chr22	2012	71140405	2012	2013
+chrX	2012	71142424	2012	2013
+chrY	2012	71144443	2012	2013
+chrM	2012	71146462	2012	2013

bio/paraphase/wrapper.py

@@ -0,0 +1,163 @@
+import glob
+import sys
+from snakemake.shell import shell
+from tempfile import TemporaryDirectory
+
+log = snakemake.log_fmt_shell(stdout=True, stderr=True)
+extra = snakemake.params.get("extra", "")
+
+
+with TemporaryDirectory() as tmpdirname:
+
+	### Run paraphase
+	shell(
+	    f"""
+	    (paraphase --bam {snakemake.input.bam} \
+	    --reference {snakemake.input.fasta} \
+	    --out {tmpdirname} \
+	    {snakemake.params.genome} \
+	    {extra}) {log}
+	    """,
+	    tmpdirname=tmpdirname,  # Pass tmpdirname to the shell environment
+	)
+	shell(
+	    """
+	    touch {snakemake.output}
+	    """
+	)
+
+	### Create a new VCF header
+
+	# Get the paths from Snakemake input and output objects
+	input_faidx = snakemake.input.faidx
+	output_vcf_header = snakemake.output.vcf_header
+
+	# Open the .fai index file and read lines
+	with open(input_faidx, "r") as fai_file:
+	    lines = fai_file.readlines()
+
+	# Open the output file and write formatted header lines
+	with open(output_vcf_header, "w") as output:
+	    for line in lines:
+		contig_id, length = line.split()[
+		    :2
+		]  # Assuming the first two elements are ID and length
+		output.write(f"##contig=<ID={contig_id},length={length}>\n")
+	output.close()
+
+	### Concatenating, reheadering, and sorting the zipped and indexed VCF files, and copy the remapped reads
+	vcf_res = glob.glob(f"{tmpdirname}/*_vcfs/*vcf")
+	if vcf_res:
+	    for vcf in vcf_res:
+		bgzip_cmd = f"bgzip -c {vcf} > {vcf}.gz"
+		shell(bgzip_cmd)
+		index_cmd = f"bcftools index {vcf}.gz"
+		shell(index_cmd)
+		print(f"Compressed and indexed: {vcf}.gz")
+
+	    params_variant_files = " ".join([f"{vcf}.gz" for vcf in vcf_res])
+	    shell(
+		f"bcftools concat -a -Oz {params_variant_files} | "
+		f"bcftools annotate --header-lines {snakemake.output.vcf_header} | "
+		f"bcftools sort -Oz -o {snakemake.output.merged_vcf}"
+	    )
+	    print(
+		f"Merged, reheadered, and sorted VCF file created: {snakemake.output.merged_vcf}"
+	    )
+
+	    # Copy out bam and bai files
+	    bam_res = glob.glob(f"{tmpdirname}/*.bam")
+	    bai_res = glob.glob(f"{tmpdirname}/*.bai")
+	    # print("BAM RES: ", bam_res, bai_res)
+	    shell(
+		f"""
+		cp -pr {' '.join(bam_res)} {snakemake.output.bam};
+		cp -pr {' '.join(bai_res)} {snakemake.output.bai}
+		"""
+	    )
+	else:
+	    print(
+		"No output VCF or BAM files were produced by paraphase, I hope this is what you were expecting, human?"
+	    )
+	    shell(f"touch {snakemake.output.merged_vcf}")
+
+
+
+
+
+
+    with TemporaryDirectory() as tmpdirname:
+
+        ### Run paraphase
+        shell(
+            f"""
+            (paraphase --bam {snakemake.input.bam} \
+            --reference {snakemake.input.fasta} \
+            --out {tmpdirname} \
+            {snakemake.params.genome} \
+            {extra}) {log}
+            """,
+            tmpdirname=tmpdirname,  # Pass tmpdirname to the shell environment
+        )
+        shell(
+            """
+            touch {snakemake.output}
+            """
+        )
+
+        ### Create a new VCF header
+
+        # Get the paths from Snakemake input and output objects
+        input_faidx = snakemake.input.faidx
+        output_vcf_header = snakemake.output.vcf_header
+
+        # Open the .fai index file and read lines
+        with open(input_faidx, "r") as fai_file:
+            lines = fai_file.readlines()
+
+        # Open the output file and write formatted header lines
+        with open(output_vcf_header, "w") as output:
+            for line in lines:
+                contig_id, length = line.split()[
+                    :2
+                ]  # Assuming the first two elements are ID and length
+                output.write(f"##contig=<ID={contig_id},length={length}>\n")
+        output.close()
+
+        ### Concatenating, reheadering, and sorting the zipped and indexed VCF files, and copy the remapped reads
+        vcf_res = glob.glob(f"{tmpdirname}/*_vcfs/*vcf")
+        if vcf_res:
+            for vcf in vcf_res:
+                bgzip_cmd = f"bgzip -c {vcf} > {vcf}.gz"
+                shell(bgzip_cmd)
+                index_cmd = f"bcftools index {vcf}.gz"
+                shell(index_cmd)
+                print(f"Compressed and indexed: {vcf}.gz")
+
+            params_variant_files = " ".join([f"{vcf}.gz" for vcf in vcf_res])
+            shell(
+                f"bcftools concat -a -Oz {params_variant_files} | "
+                f"bcftools annotate --header-lines {snakemake.output.vcf_header} | "
+                f"bcftools sort -Oz -o {snakemake.output.merged_vcf}"
+            )
+            print(
+                f"Merged, reheadered, and sorted VCF file created: {snakemake.output.merged_vcf}"
+            )
+
+            # Copy out bam and bai files
+            bam_res = glob.glob(f"{tmpdirname}/*.bam")
+            bai_res = glob.glob(f"{tmpdirname}/*.bai")
+            # print("BAM RES: ", bam_res, bai_res)
+            shell(
+                f"""
+                cp -pr {' '.join(bam_res)} {snakemake.output.bam};
+                cp -pr {' '.join(bai_res)} {snakemake.output.bai}
+                """
+            )
+        else:
+            print(
+                "No output VCF or BAM files were produced by paraphase, I hope this is what you were expecting, human?"
+            )
+            shell(f"touch {snakemake.output.merged_vcf}")
+
+