add syndna command (#15)

* add syndna command

* fix error

* adding tests

* syndna CONDA_ENVIRONMENT

Author: antgonza

.github/workflows/qiita-plugin-ci.yml

@@ -145,7 +145,7 @@ jobs:
   lint:
     runs-on: ubuntu-latest
     steps:
-      - name: flake8
+      - name: ruff
         uses: actions/setup-python@v2
         with:
           python-version: 3.9
@@ -155,5 +155,5 @@ jobs:
         uses: actions/checkout@v2
       - name: lint
         run: |
-          pip install -q flake8
-          flake8 .
+          pip install -q ruff
+          ruff check .

README.rst

@@ -8,11 +8,11 @@ Qiita plugin to process PacBio reads; it currently provides 2 commands for Qiita
 * **PacBio processing**: which goes from step 1 to 7 in the image below. The expected output
   is a main folder with folders per-sample and folders for each of the different outputs, as follows:
 
-  * **MAG** folder: all Metagenome-Assembled Genome (MAG) generatedfor that sample
+  * **MAG** folder: all Metagenome-Assembled Genome (MAG) generated for that sample
   * **LCG** folder: all Long-Circular Genome (LCG) generated for that sample that are over 512kb in size - approximate 515,000 bases (half a million)
   * **small_LCG** folder: all Long-Circular Genome (LCG) generated for that sample that are under 512kb in size
   * **[sample-name].fna.gz**: the no LCG reads used for MAG generation
-  * **[sample-name].checkm.txt.gz**: MAG quanlity information from CheckM v1.2.3
+  * **[sample-name].checkm.txt.gz**: MAG quality information from CheckM v1.2.3
 
 
 .. image:: images/PacBioProcessing.png

pyproject.toml

@@ -36,7 +36,7 @@ dependencies = [
     "click>=3.3",
     "pandas",
     'requests',
-    'flake8',
+    'ruff',
     'coverage',
     'pytest-cov',
     'numpy',
@@ -46,6 +46,8 @@ dependencies = [
     "qiita-files@https://github.com/qiita-spots/qiita-files/archive/master.zip",
     "qiita_client@https://github.com/qiita-spots/qiita_client/archive/master.zip",
     "woltka@git+https://github.com/qiyunzhu/woltka.git#egg=woltka",
+    "pysyndna@git+https://github.com/biocore/pysyndna.git#egg=pysyndna",
+    'micov',
 ]
 
 [project.scripts]

qp_pacbio/__init__.py

@@ -5,10 +5,14 @@
 #
 # The full license is in the file LICENSE, distributed with this software.
 # -----------------------------------------------------------------------------
-from qiita_client import QiitaPlugin, QiitaCommand
-from .qp_pacbio import pacbio_processing, minimap2_processing
-from .util import plugin_details
+from qiita_client import QiitaCommand, QiitaPlugin
 
+from .qp_pacbio import (
+    minimap2_processing,
+    pacbio_processing,
+    syndna_processing,
+)
+from .util import plugin_details
 
 plugin = QiitaPlugin(**plugin_details)
 
@@ -60,3 +64,29 @@
 )
 
 plugin.register_command(pacbio_processing_cmd)
+
+#
+# syndna filtering
+#
+
+req_params = {
+    "artifact": ("artifact", ["per_sample_FASTQ"]),
+}
+opt_params = {"min_sample_counts": ("integer", "1")}
+outputs = {
+    "SynDNA hits": "BIOM",
+    "reads without SynDNA": "per_sample_FASTQ",
+}
+dflt_param_set = {
+    "SynDNA": {"min_sample_counts": 1},
+}
+pacbio_processing_cmd = QiitaCommand(
+    "Remove SynDNA plasmid, insert, & GCF_000184185 reads (minimap2)",
+    "Remove SynDNA reads using minimap2 and woltka",
+    syndna_processing,
+    req_params,
+    opt_params,
+    outputs,
+    dflt_param_set,
+)
+plugin.register_command(pacbio_processing_cmd)

qp_pacbio/data/resources.yaml

@@ -58,3 +58,15 @@ Woltka v0.1.7, minimap2:
     nprocs: 16
     wall_time_limit: 1-00:00:00
     mem_in_gb: 120
+Remove SynDNA plasmid, insert, & GCF_000184185 reads (minimap2):
+  syndna:
+    node_count: 1
+    nprocs: 16
+    wall_time_limit: 10:00:00
+    mem_in_gb: 60
+    max_tasks: 16
+  finish:
+    node_count: 1
+    nprocs: 16
+    wall_time_limit: 1-00:00:00
+    mem_in_gb: 120

qp_pacbio/data/templates/syndna.sbatch

@@ -23,25 +23,28 @@ sample_name=`echo $input | awk '{print $1}'`
 filename=`echo $input | awk '{print $2}'`
 fn=`basename ${filename}`
 
-mkdir -p ${out_folder}/clean/
+mkdir -p ${out_folder}/filtered/
+
+sn_folder=${out_folder}/bioms/${sample_name}
+mkdir -p ${sn_folder}
+
+coverm contig --single $filename --reference ${db_folder}/All_synDNA_inserts.fasta --mapper minimap2-hifi \
+    --min-read-percent-identity 0.95 --min-read-aligned-percent 0.0 -m mean count --threads {{nprocs}} \
+    --output-file ${sn_folder}/${sample_name}.txt
+cat ${sn_folder}/${sample_name}_insert_counts.txt | sed 's/Contig/\#OTU ID/' | \
+    sed 's/ Read Count//' > ${sn_folder}/${sample_name}.tsv
+biom convert -i ${sn_folder}/${sample_name}.txt -o ${sn_folder}/${sample_name}.biom --to-hdf5
 
 # removing AllsynDNA_plasmids_FASTA_ReIndexed_FINAL.fasta not coverm
 minimap2 -x map-hifi -t {{nprocs}} -a --MD --eqx -o ${out_folder}/${sample_name}_plasmid.sam ${db_folder}/AllsynDNA_plasmids_FASTA_ReIndexed_FINAL.fasta $filename
 samtools view -F 4 -@ {{nprocs}} ${out_folder}/${sample_name}_plasmid.sam | awk '{print $1}' | sort -u > ${out_folder}/${sample_name}_plasmid_mapped.txt
 seqkit grep -v -f ${out_folder}/${sample_name}_plasmid_mapped.txt $filename > ${out_folder}/${sample_name}_no_plasmid.fastq
 
-# removing All_synDNA_inserts.fasta use coverm
-minimap2 -x map-hifi -t {{nprocs}} -a --MD --eqx -o ${out_folder}/${sample_name}_inserts.sam ${db_folder}/All_synDNA_inserts.fasta ${out_folder}/${sample_name}_no_plasmid.fastq
-samtools view -bS -@ {{int(nprocs/2)}} ${out_folder}/${sample_name}_no_plasmid.fastq | samtools sort -@ {{int(nprocs/2)}} -O bam -o ${out_folder}/${sample_name}_inserts_sorted.sam
-coverm filter --bam-files ${out_folder}/${sample_name}_inserts_sorted.sam --min-read-percent-identity 99.9 --min-read-aligned-percent 95 --threads {{nprocs}} -o ${out_folder}/${sample_name}_no_inserts.bam
-samtools view -O SAM -o ${out_folder}/${sample_name}_no_inserts_sorted.sam ${out_folder}/${sample_name}_no_inserts.bam
-awk '{print $1}' ${out_folder}/${sample_name}_no_inserts_sorted.sam > ${out_folder}/${sample_name}_reads_filtered.txt
-seqkit grep -v -f ${out_folder}/${sample_name}_reads_filtered.txt ${out_folder}/${sample_name}_no_plasmid.fastq > ${out_folder}/${sample_name}_no_plasmid_no_inserts.fastq
-
 # removing GCF_000184185.1_ASM18418v1_genomic_chroso.fna use coverm
 minimap2 -x map-hifi -t {{nprocs}} -a --MD --eqx -o ${out_folder}/${sample_name}_GCF_000184185.sam ${db_folder}/GCF_000184185.1_ASM18418v1_genomic_chroso.fna ${out_folder}/${sample_name}_no_plasmid_no_inserts.fastq
-samtools view -bS -@ {{int(nprocs/2)}} ${out_folder}/${sample_name}_no_plasmid_no_inserts.fastq | samtools sort -@ {{int(nprocs/2)}} -O bam -o ${out_folder}/${sample_name}_GCF_000184185_sorted.sam
+samtools view -bS -@ {{ nprocs/2 | int }} ${out_folder}/${sample_name}_no_plasmid_no_inserts.fastq | samtools sort -@ {{ nprocs/2 | int }} -O bam -o ${out_folder}/${sample_name}_GCF_000184185_sorted.sam
 coverm filter --bam-files ${out_folder}/${sample_name}_GCF_000184185_sorted.sam --min-read-percent-identity 99.9 --min-read-aligned-percent 95 --threads {{nprocs}} -o ${out_folder}/${sample_name}_GCF_000184185.bam
 samtools view -O SAM -o ${out_folder}/${sample_name}_no_GCF_000184185_sorted.sam ${out_folder}/${sample_name}_no_inserts.bam
 awk '{print $1}' ${out_folder}/${sample_name}_no_GCF_000184185_sorted.sam > ${out_folder}/${sample_name}_GCF_000184185_reads_filtered.txt
-seqkit grep -v -f ${out_folder}/${sample_name}_GCF_000184185_reads_filtered.txt ${out_folder}/${sample_name}_GCF_000184185.fastq | gz > ${out_folder}/clean/${fn}
+seqkit grep -v -f ${out_folder}/${sample_name}_GCF_000184185_reads_filtered.txt ${out_folder}/${sample_name}_GCF_000184185.fastq | gz > ${out_folder}/filtered/${fn}
+awk 'BEGIN {FS=OFS="\t"}; {print $1,$3}'

qp_pacbio/data/templates/syndna_finish.sbatch

@@ -0,0 +1,24 @@
+#!/bin/bash
+#SBATCH -J {{job_name}}
+#SBATCH -p qiita
+#SBATCH -N {{node_count}}
+#SBATCH -n {{nprocs}}
+#SBATCH --time {{wall_time_limit}}
+#SBATCH --mem {{mem_in_gb}}G
+#SBATCH -o {{output}}/merge/logs/%x-%A_%a.out
+#SBATCH -e {{output}}/merge/logs/%x-%A_%a.err
+
+source ~/.bashrc
+set -e
+{{conda_environment}}
+cd {{output}}/
+
+biom_merge_pacbio --base {{output}} --type syndna
+
+# find {{output}}/coverages/ -iname "*.cov" > {{output}}/cov_files.txt
+# micov consolidate --paths {{output}}/cov_files.txt --lengths ${len_map} --output {{output}}/coverages.tgz
+
+# cd alignment
+# tar -cvf ../alignment.tar *.sam.xz
+
+finish_qp_pacbio {{url}} {{qjid}} {{output}}

qp_pacbio/data/templates/woltka_minimap2_merge.sbatch

@@ -39,7 +39,7 @@ for f in `ls bioms/*/per-gene.biom`; do
 done | parallel --halt now,fail=1 -j {{nprocs}}
 wait
 
-biom_merge_pacbio --base {{output}}
+biom_merge_pacbio --base {{output}} --type woltka
 
 find {{output}}/coverages/ -iname "*.cov" > {{output}}/cov_files.txt
 micov consolidate --paths {{output}}/cov_files.txt --lengths ${len_map} --output {{output}}/coverages.tgz

qp_pacbio/qp_pacbio.py

@@ -11,8 +11,13 @@
 from shutil import copy2
 from subprocess import run
 
+import pandas as pd
 import yaml
+from biom import load_table
 from jinja2 import Environment
+from pysyndna import (
+    fit_linear_regression_models_for_qiita,
+)
 from qiita_client import ArtifactInfo
 
 from .util import KISSLoader, find_base_path
@@ -424,3 +429,147 @@ def generate_minimap2_processing(qclient, job_id, out_dir, parameters, url):
     minimap2_merge_script = _write_slurm(f"{out_dir}/merge", m2mt, **params)
 
     return minimap2_script, minimap2_merge_script
+
+
+def syndna_processing(qclient, job_id, parameters, out_dir):
+    """generates output for syndna processing.
+
+    Parameters
+    ----------
+    qclient : tgp.qiita_client.QiitaClient
+        Qiita server client.
+    job_id : str
+        Job id.
+    parameters : dict
+        Parameters for this job.
+    out_dir : str
+        Output directory.
+
+    Returns
+    -------
+    bool, list, str
+        Results tuple for Qiita.
+    """
+    qclient.update_job_step(job_id, "Commands finished")
+
+    errors = []
+    ainfo = []
+    fp_biom = f"{out_dir}/syndna.biom"
+    # do we need to stor alignments?
+    # fp_alng = f'{out_dir}/sams/final/alignment.tar'
+
+    if exists(fp_biom):  # and exists(fp_alng):
+        # if we got to this point a preparation file should exist in
+        # the output folder
+        prep = pd.read_csv(f"{out_dir}/prep_info.tsv", index_col=None, sep="\t")
+        output = fit_linear_regression_models_for_qiita(
+            prep, load_table(fp_biom), int(parameters["min_sample_counts"])
+        )
+        # saving results to disk
+        lin_regress_results_fp = f"{out_dir}/lin_regress_by_sample_id.yaml"
+        fit_syndna_models_log_fp = f"{out_dir}/fit_syndna_models_log.txt"
+        with open(lin_regress_results_fp, "w") as fp:
+            fp.write(output["lin_regress_by_sample_id"])
+        with open(fit_syndna_models_log_fp, "w") as fp:
+            fp.write(output["fit_syndna_models_log"])
+        ainfo = [
+            ArtifactInfo(
+                "SynDNA hits",
+                "BIOM",
+                [
+                    (fp_biom, "biom"),
+                    # rm if fp_alng is not needed
+                    # (fp_alng, "log"),
+                    (lin_regress_results_fp, "log"),
+                    (fit_syndna_models_log_fp, "log"),
+                ],
+            )
+        ]
+    else:
+        ainfo = []
+        errors.append(
+            'Missing files from the "SynDNA hits"; please '
+            "contact [email protected] for more information"
+        )
+
+    fp_seqs = f"{out_dir}/filtered"
+    reads = []
+    for f in glob(f"{fp_seqs}/*.fastq.gz"):
+        reads.append((f, "raw_forward_seqs"))
+
+    if not errors:
+        ainfo.append(ArtifactInfo("reads without SynDNA", "per_sample_FASTQ", reads))
+    else:
+        return False, ainfo, "\n".join(errors)
+
+    return True, ainfo, ""
+
+
+def generate_syndna_processing(qclient, job_id, out_dir, parameters, url):
+    """generates slurm scripts for syndna processing.
+
+    Parameters
+    ----------
+    qclient : tgp.qiita_client.QiitaClient
+        Qiita server client.
+    job_id : str
+        Job id.
+    out_dir : str
+        Output directory.
+    parameters : dict
+        Parameters for this job.
+    url : str
+        URL to send the respose, finish the job.
+
+    Returns
+    -------
+    str, str
+        Returns the two filepaths of the slurm scripts
+    """
+    resources = RESOURCES[
+        "Remove SynDNA plasmid, insert, & GCF_000184185 reads (minimap2)"
+    ]
+    main_parameters = {
+        "conda_environment": CONDA_ENVIRONMENT,
+        "output": out_dir,
+        "qjid": job_id,
+    }
+
+    qclient.update_job_step(
+        job_id, "Step 1 of 4: Collecting info and generating submission"
+    )
+
+    artifact_id = parameters["artifact"]
+
+    njobs = generate_sample_list(qclient, artifact_id, out_dir)
+
+    qclient.update_job_step(
+        job_id,
+        "Step 2 of 4: Creating submission templates",
+    )
+
+    m2t = JGT("syndna.sbatch")
+    step_resources = resources["syndna"]
+    params = main_parameters | {
+        "job_name": f"sd_{job_id}",
+        "node_count": step_resources["node_count"],
+        "nprocs": step_resources["nprocs"],
+        "wall_time_limit": step_resources["wall_time_limit"],
+        "mem_in_gb": step_resources["mem_in_gb"],
+        "array_params": f"1-{njobs}%{step_resources['max_tasks']}",
+    }
+    minimap2_script = _write_slurm(f"{out_dir}/minimap2", m2t, **params)
+
+    m2mt = JGT("syndna_finish.sbatch")
+    step_resources = resources["finish"]
+    params = main_parameters | {
+        "job_name": f"me_{job_id}",
+        "node_count": step_resources["node_count"],
+        "nprocs": step_resources["nprocs"],
+        "wall_time_limit": step_resources["wall_time_limit"],
+        "mem_in_gb": step_resources["mem_in_gb"],
+        "url": url,
+    }
+    minimap2_merge_script = _write_slurm(f"{out_dir}/merge", m2mt, **params)
+
+    return minimap2_script, minimap2_merge_script