analysis: FAA files for metaP analysis (all, cluster representatives) (issue #97)

workflow/envs/biopython.yaml

+../../workflow_zymo/envs/biopython.yaml
\ No newline at end of file

workflow/rules/analysis.smk

@@ -363,6 +363,20 @@ rule analysis_mmseqs2:
         "mmseqs createtsv {params.db} {params.db} {params.cl} {output.tsv} && "
         "date) &> {log}"
 
+# FAA file of cluster representatives
+rule analysis_mmseqs2_cluster_faa:
+    input:
+        faa=os.path.join(RESULTS_DIR, "analysis/mmseqs2/proteins.faa"),
+        tsv=os.path.join(RESULTS_DIR, "analysis/mmseqs2/clusters.tsv")
+    output:
+        os.path.join(RESULTS_DIR, "analysis/mmseqs2/clusters.faa")
+    conda:
+        os.path.join(ENV_DIR, "biopython.yaml")
+    message:
+        "Analysis: MMSeqs2 clustering FAA: {input}"
+    script:
+        os.path.join(SRC_DIR, "mmseqs2_cluster_faa.py")
+
 ##################################################
 # COVERAGE
 

workflow/scripts/mmseqs2_cluster_faa.py

+#!/usr/bin/env python
+
+# FAA file containing only cluster representatives
+
+import pandas
+from Bio import SeqIO
+
+# protein IDs = cluster representatives (1st column)
+prot_ids = set(pandas.read_csv(snakemake.input.tsv, sep="\t", header=None).iloc[:,0])
+# filter input FAA by IDs
+with open(snakemake.input.faa, "r") as ifile, open(snakemake.output[0], "w") as ofile:
+    for record in SeqIO.parse(ifile, "fasta"):
+        if record.id in prot_ids:
+            SeqIO.write(record, ofile, "fasta")

workflow/steps/analysis.smk

@@ -49,6 +49,10 @@ rule ANALYSIS:
         ) if "diamond" in config["steps_analysis"] else [],
         # mmseqs2
         [os.path.join(RESULTS_DIR, "analysis/mmseqs2/clusters.tsv")] if "mmseqs2" in config["steps_analysis"] else [],
+        [ # for metaP analysis
+            os.path.join(RESULTS_DIR, "analysis/mmseqs2/proteins.faa"),
+            os.path.join(RESULTS_DIR, "analysis/mmseqs2/clusters.faa")
+        ] if "mmseqs2" in config["steps_analysis"] and "metap" in config["data"] else [],
         # coverage
         expand(
             os.path.join(RESULTS_DIR, "mapping/{mtype}/{rtype_tool}.cov.pergene"),