run_ipsn.snake

#!/usr/bin/env python

# Pipeline to extract proteins from assembly with prokka and obtain functional annotations with interproscan
import os
import re
import glob
import sys

# Make directories if they don't already exist
os.system('mkdir -pv genomes prokka interproscan/kegg interproscan/metacyc interproscan/reactome/ interproscan/go_bp interproscan/go_mf interproscan/go_cc interproscan/json_output interproscan/dirty')

# Find list of input gzipped fna files:
genome_files = glob.glob('genomes/*fna*gz')

# Get list of final output files
clean_files = []

for file in genome_files:
    file = os.path.basename(file)
    file = file.split('.fna')[0]
    reactome_file = ("interproscan/reactome/%s_reactome.txt" %(file))
    kegg_file = ("interproscan/kegg/%s_kegg.txt" %(file))
    metacyc_file = ("interproscan/metacyc/%s_metacyc.txt" %(file))
    bp_file = ("interproscan/go_bp/%s_bp.txt" %(file))
    mf_file = ("interproscan/go_mf/%s_mf.txt" %(file))
    cc_file= ("interproscan/go_cc/%s_cc.txt" %(file))
    clean_files.append(reactome_file)
    clean_files.append(kegg_file)
    clean_files.append(metacyc_file)
    clean_files.append(bp_file)
    clean_files.append(mf_file)
    clean_files.append(cc_file)
#print(clean_files)

rule all:
    input: [file for file in clean_files]

rule clean_reactome:
    input: "interproscan/dirty/{sample}_reactome.dirty"
    output: "interproscan/reactome/{sample}_reactome.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule clean_kegg:
    input: "interproscan/dirty/{sample}_kegg.dirty"
    output: "interproscan/kegg/{sample}_kegg.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule clean_metacyc:
    input: "interproscan/dirty/{sample}_metacyc.dirty"
    output: "interproscan/metacyc/{sample}_metacyc.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule clean_bp:
    input: "interproscan/dirty/{sample}_bp.dirty"
    output: "interproscan/go_bp/{sample}_bp.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule clean_mf:
    input: "interproscan/dirty/{sample}_mf.dirty"
    output: "interproscan/go_mf/{sample}_mf.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule clean_cc:
    input: "interproscan/dirty/{sample}_cc.dirty"
    output: "interproscan/go_cc/{sample}_cc.txt"
    priority: 50
    shell: "bin/clean_pathways.py -i {input} -o {output}"

rule extract_reactome:
    input: json_in="interproscan/json/{sample}.json"
    output: reactome="interproscan/dirty/{sample}_reactome.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 1 'Reactome'  | grep 'name' > {output.reactome}"

rule extract_kegg:
    input: json_in="interproscan/json/{sample}.json"
    output: kegg="interproscan/dirty/{sample}_kegg.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 1 'KEGG'  | grep 'name' > {output.kegg}"

rule extract_metacyc:
    input: json_in="interproscan/json/{sample}.json"
    output: metacyc="interproscan/dirty/{sample}_metacyc.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 1 'MetaCyc'  | grep 'name' > {output.metacyc}"

rule extract_bp:
    input: json_in="interproscan/json/{sample}.json"
    output: bp="interproscan/dirty/{sample}_bp.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 2 'BIOLOGICAL_PROCESS' | grep 'name' > {output.bp}"

rule extract_mf:
    input: json_in="interproscan/json/{sample}.json"
    output: mf="interproscan/dirty/{sample}_mf.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 2 'MOLECULAR_FUNCTION' | grep 'name' > {output.mf}"

rule extract_cc:
    input: json_in="interproscan/json/{sample}.json"
    output: cc="interproscan/dirty/{sample}_cc.dirty"
    priority: 40
    shell: "cat {input.json_in} | grep -B 2 'CELLULAR_COMPONENT' | grep 'name' > {output.cc}"

rule interproscan:
    input: protein_in="prokka/{sample}/{sample}.faa"
    output: ipscn_out="interproscan/json/{sample}.json" 
    priority: 30
    shell: "singularity exec bin/sepathway.sif interproscan.sh -f json --goterms --iprlookup --pathways -T temp -i {input.protein_in} -o {output.ipscn_out}"

rule prokka_annotate:
    input:  genome_in="genomes/{sample}.fna"
    output: "prokka/{sample}/{sample}.faa"
    priority: 12
    shell: "SINGULARITYENV_LC_ALL=C singularity exec --cleanenv bin/prokka.sif prokka {input.genome_in} --outdir prokka/{wildcards.sample} --prefix {wildcards.sample} --quiet --force ; gzip {input.genome_in}"

rule gunzip_genomes:
    input: i1="genomes/{sample}.fna.gz" 
    output: o1="genomes/{sample}.fna"
    priority: 10
    shell: "gunzip {input.i1}"