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,15 @@
+name: poraphrase
+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,143 @@
+import glob
+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")
+
+    ### 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}")