Snakefile.library

# Obtain run_ids from DDA or DIA-ps mzXML files in data folder
run_ids, = glob_wildcards("data_dda/{run}.mzXML")

rule all:
    input:
        global_target_pqp =  expand("results/library/{run}_global_pqp.tsv", run=run_ids),
        pqp = "data_library/library.pqp"

rule decoydb:
    input:
        "data_library/library.fasta"
    output:
        "results/library/library.fasta"
    singularity:
        "docker://openswath/develop:latest"
    shell:
        "DecoyDatabase -in {input} -out {output}"

rule digestdb:
    input:
        rules.decoydb.output
    output:
        "results/library/peptides.fasta"
    singularity:
        "docker://openswath/develop:latest"
    shell:
        "Digestor -in {input} -out {output} -missed_cleavages 1 -min_length 7 -max_length 50 -enzyme Trypsin"

rule peptide_index:
    input:
        rules.digestdb.output
    output:
        "results/library/peptides.pkl"
    singularity:
        "docker://grosenberger/easypqp:latest"
    shell:
        "easypqp index --fasta {input} --pepidx {output}"

rule msfragger_index:
    input:
        fasta = rules.decoydb.output,
        mzxml = "params/msfragger.mzXML"
    output:
        index = "results/library/library.fasta.1.pepindex",
        tmp = temp("params/msfragger.tsv")
    params:
        fragger = "params/fragger_closed.params"
    singularity:
        "docker://grosenberger/msfragger:20190222"
    threads: 4
    resources:
        mem_mb=lambda wildcards, attempt: (attempt * 32768) - 2048
    shell:
        "java -Xmx{resources.mem_mb}M -jar /MSFragger.jar {params.fragger} {input.mzxml}"

rule msfragger_search:
    input:
        fasta = rules.decoydb.output,
        index = rules.msfragger_index.output.index,
        mzxml = "data_dda/{run}.mzXML"
    output:
        cache = temp("results/library/{run}.mzXML"),
        fragger = "results/library/{run}.tsv"
    params:
        fragger = "params/fragger_closed.params"
    singularity:
        "docker://grosenberger/msfragger:20190222"
    threads: 4
    resources:
        mem_mb=lambda wildcards, attempt: (attempt * 32768) - 2048
    shell:
        "ln {input.mzxml} {output.cache} && "
        "java -Xmx{resources.mem_mb}M -jar /MSFragger.jar {params.fragger} {output.cache}"

rule msfragger_convert:
    input:
        peptides = rules.peptide_index.output,
        mzxml = "data_dda/{run}.mzXML",
        fragger = rules.msfragger_search.output.fragger
    output:
        pyprophet = "results/library/{run}_pyprophet.tsv",
        subsampled = "results/library/{run}_subsampled.tsv",
        peakpkl = "results/library/{run}.peakpkl"
    params:
        unimod = "params/unimod_phospho.xml",
	subsample_fraction = 1.0
    singularity:
        "docker://grosenberger/easypqp:latest"
    shell:
        "easypqp convert --unimod {params.unimod} --fragger {input.fragger} --pepidx {input.peptides} --psms {output.pyprophet} --subpsms {output.subsampled} --mzxml {input.mzxml} --peaks {output.peakpkl} --subsample_fraction {params.subsample_fraction}"

rule pyprophet_learn:
    input:
        expand("results/library/{run}_subsampled.tsv", run=run_ids),
    output:
        "results/library/pyprophet_learn_ms2_model.bin"
    params:
        merged = "results/library/pyprophet_learn.tsv",
        subsample_factor = 1
    singularity:
        "docker://pyprophet/master:latest"
    threads: 4
    shell:
        "awk 'BEGIN {{ FS=\"\t\"; OFS=\"\t\" }} FNR>1 || NR==1 {{print $0}}' {input} > {params.merged} && "
        "pyprophet score --in {params.merged} --threads={threads} --classifier=XGBoost --xeval_num_iter=3 --ss_initial_fdr=0.1 --ss_iteration_fdr=0.05"

rule pyprophet_score:
    input:
        pyprophet =rules.msfragger_convert.output.pyprophet,
        model = rules.pyprophet_learn.output
    output:
        "results/library/{run}_pyprophet_scored.tsv"
    singularity:
        "docker://pyprophet/master:latest"
    shell:
        "pyprophet score --in {input.pyprophet} --classifier=XGBoost --apply_weights={input.model}"

rule easypqp:
    input:
        psms = expand("results/library/{run}_pyprophet_scored.tsv", run=run_ids),
        peakpkl = expand("results/library/{run}.peakpkl", run=run_ids),
    output:
        peptide_plot = "results/library/pyprophet_peptides.pdf",
        protein_plot = "results/library/pyprophet_protein.pdf",
    singularity:
        "docker://grosenberger/easypqp:latest"
    params:
        psm_fdr_threshold = 0.01,
        peptide_fdr_threshold = 0.01,
        protein_fdr_threshold = 0.01
    shell:
        "easypqp library --psm_fdr_threshold={params.psm_fdr_threshold} --peptide_fdr_threshold={params.peptide_fdr_threshold} --protein_fdr_threshold={params.protein_fdr_threshold} --peptide_plot={output.peptide_plot} --protein_plot={output.protein_plot} {input.psms} {input.peakpkl}"

rule global_target_pqp:
    input:
        iRT = rules.easypqp.output
    params:
        unimod = "params/unimod_phospho.xml",
        peaks = temp("results/library/{run}_global_peaks.tsv")
    output:
        temp("results/library/{run}_global_pqp.tsv")
    singularity:
        "docker://openswath/develop:latest"
    shell:
        #"OpenSwathAssayGenerator -in {params.peaks} -out {output}"
        "OpenSwathAssayGenerator -in {params.peaks} -out {output} -enable_ipf -unimod_file {params.unimod} -disable_identification_ms2_precursors -disable_identification_specific_losses"

rule global_combined_pqp:
    input:
        expand("results/library/{run}_global_pqp.tsv", run=run_ids)
    output:
        temp("results/library/combined_global_pqp.tsv")
    shell:
        "awk 'BEGIN {{ FS=\"\t\"; OFS=\"\t\" }} FNR>1 || NR==1 {{print $1,$2,$3,$4,$5,$6,$7,$8,$15,$18,$19,$20,$25,$26,$27,$28,$29}}' {input} > {output}"

rule global_combined_decoy_pqp:
    input:
        rules.global_combined_pqp.output
    output:
        "data_library/library.pqp"
    singularity:
        "docker://openswath/develop:latest"
    shell:
        "cache=${{TMPDIR-/tmp/}} && "
        "OpenSwathDecoyGenerator -in {input} -out $cache/library.pqp && mv $cache/library.pqp {output}"