neoflow_msms.nf

#!/usr/bin/env nextflow


/*
 * Default pipeline parameters. They can be overriden on the command line eg.
 * given `params.ms` specify on the run command line `--ms ms_data_folder`.
 */

params.help = false
params.pv_enzyme = 1
params.pv_c      = 2
params.pv_tol    = 10
params.pv_tolu   = "ppm"
params.pv_itol   = 0.5 // Da
params.pv_fixmod = 6
params.pv_varmod = 107

params.search_engine = "comet"
params.rt_validation = false
params.prefix = "pga"
params.out_dir = "./"



/* Prints help when asked for and exits */

def helpMessage() {
    log.info"""
    =========================================
    neoflow => Variant peptide identification
    =========================================
    Usage:
    nextflow run neoflow-msms.nf
    MS/MS searching arguments:
      --db                        The customized protein database (target + decoy sequences) in FASTA format which is generated by neoflow_db.nf
      --ms                        MS/MS data in MGF format
      --msms_para_file            Parameter file for MS/MS searching
      --out_dir                   Output folder, default is "./"
      --prefix                    The prefix of output files
      --search_engine             The search engine used for MS/MS searching, comet=Comet, msgf=MS-GF+ or xtandem=X!Tandem

    PepQuery arguments:
      --pv_enzyme                 Enzyme used for protein digestion. 0:Non enzyme, 1:Trypsin (default), 2:Trypsin (no P rule), 3:Arg-C, 4:Arg-C (no P rule), 5:Arg-N, 6:Glu-C, 7:Lys-C
      --pv_c                      The max missed cleavages, default is 2
      --pv_tol                    Precursor ion m/z tolerance, default is 10
      --pv_tolu                   The unit of --tol, ppm or Da. Default is ppm
      --pv_itol                   The error window for fragment ion, default is 0.5
      --pv_fixmod                 Fixed modification. The format is like : 1,2,3. Different modification is represented by different number
      --pv_varmod                 Variable modification. The format is the same with --fixMod;
      --pv_refdb                  Reference protein database

    AutoRT parameters:
      --rt_validation             Perform RT based validation                     

	  --help                      Print help message
    """.stripIndent()
}


// Show help emssage
if (params.help){
    helpMessage()
    exit 0
}

// input
//msms_dir = file(params.msms_dir)
msms_para_file = file(params.msms_para_file)
search_engine = params.search_engine
rt_validation = params.rt_validation
search_db = file(params.db)
out_dir = file(params.out_dir)
// a folder which contains multiple MGF files or a single MGF file
ms_data = file(params.ms)
out_prefix = params.prefix

// PepQuery parameters
pv_enzyme = params.pv_enzyme
pv_c      = params.pv_c
pv_tol    = params.pv_tol
pv_tolu   = params.pv_tolu
pv_itol   = params.pv_itol
pv_fixmod = params.pv_fixmod
pv_varmod = params.pv_varmod
pv_refdb = file(params.pv_refdb)

if(!out_dir.isDirectory()){
    out_dir_result = out_dir.mkdirs()
    println out_dir_result ? "Create folder: $out_dir!" : "Cannot create directory: $myDir!"
}



if(ms_data.isFile()){
    println "Process single MS/MS file."
    ms_data_file = file(params.ms)
}else{
    println "Process multiple MS/MS files."
    ms_data_file = Channel.fromPath("${params.ms}/*.mgf")
}


process msms_searching{

    tag "${ms_file}"

    container "proteomics/neoflow:latest"
    
    cpus "$params.cpu"

    publishDir "${out_dir}/msms_searching/", mode: "copy", overwrite: true

    input:
    file(ms_file) from ms_data_file
    file(msms_para_file)
    file(search_db)

    output:
    file("${res_file}") into psm_raw_files
    

    script:

    if (search_engine == "msgf") {
        res_file = "${ms_file.baseName}.mzid"
        """
        #!/bin/sh
        java -Xmx10g -jar /opt/MSGFPlus.jar \
            -s ${ms_file} \
            -d ${search_db} \
            -conf ${msms_para_file} \
            -tda 0 \
            -o ${ms_file.baseName}.mzid
        """
    }else if(search_engine == "comet") {
        res_file = "${ms_file.baseName}_rawResults.txt"
        """
        #!/bin/sh
        /opt/comet.2018014.linux.exe -P${msms_para_file} -N${ms_file.baseName}_rawResults -D${search_db} ${ms_file}
        sed -i '1d' ${ms_file.baseName}_rawResults.txt
        sed -i '1 s/\$/\tna/' ${ms_file.baseName}_rawResults.txt
        """
    }else if(search_engine == "xtandem"){

        ms_file_name = "${ms_file.baseName}"
        xml_input = "${ms_file.baseName}_input.xml"
        res_file = "${ms_file_name}.mzid"
        """
        #!/bin/sh
        python3 ${baseDir}/bin/generate_xtandem_para_xml.py ${msms_para_file} ${ms_file} ${search_db} ${ms_file_name} ${xml_input}
        
        ## users must provide the main search parameter file for X!Tandem search
        /opt/tandem-linux-17-02-01-4/bin/tandem.exe ${xml_input}

        ## convert xml to mzid
        java -Xmx10g -jar /opt/mzidlib-1.7/mzidlib-1.7.jar Tandem2mzid \
            ${ms_file_name}.xml \
            ${ms_file_name}.mzid \
            -outputFragmentation false \
            -decoyRegex XXX_ \
            -databaseFileFormatID MS:1001348 \
            -massSpecFileFormatID MS:1001062 \
            -idsStartAtZero false \
            -compress false \
            -proteinCodeRegex "\\S+"
        """
    }

}


process calculate_fdr{

    tag "calculate_fdr"

    container "proteomics/pga:latest"
    
    cpus "$params.cpu"

    publishDir "${out_dir}/fdr_estimation/", mode: "copy", overwrite: true

    input:
    file(psm_raw_file) from psm_raw_files.collect()
    file(search_db)

    output:
    file("*_level") into (pga_result_folder,pga_result_folder_for_autort)
    file("*-rawPSMs.txt") into raw_psm_file

    script:
    """
    mkdir -p peptide_level
    mkdir -p psm_level
    mkdir -p peptide_level/global_fdr
    mkdir -p psm_level/global_fdr
    Rscript ${baseDir}/bin/fdr_calc.R ./ ${search_db} ${out_prefix} ${search_engine} ./
    """
    
}


process prepare_pepquery_input{

    tag "prepare_pepquery_input"

    container "proteomics/pga:latest"
    
    cpus "$params.cpu"

    input:
    file(pga_result_folder)

    output:
    //
    file("novel_peptides_psm.tsv") into novel_psm_tsv
    file("novel_peptides.tsv") into novel_peptide_tsv

    script:
    """
    #!/usr/bin/env /usr/local/bin/Rscript
    library(dplyr)
    library(readr)
    library(stringr)

    is_novel_peptides=function(protein_ids){
        ids <- str_split(protein_ids,pattern=";") %>% unlist()
        is_novel <- all(str_detect(ids,pattern="^VAR"))
        return(is_novel)
    }
    
    psms <- read_tsv("peptide_level/global_fdr/pga-peptideSummary.txt")


    novel_psms <- psms %>% filter(Qvalue<=0.01) %>% 
        mutate(is_novel=sapply(protein,is_novel_peptides)) %>%
        filter(is_novel==TRUE) %>%
        filter(isdecoy==FALSE)

    ## output
    out_novel_psm_file <- "novel_peptides_psm.tsv"
    novel_psms %>% write_tsv(out_novel_psm_file)

    novel_psms %>% select(peptide) %>% distinct() %>% write_tsv("novel_peptides.tsv",col_names=FALSE)
    """
}


process run_pepquery{

    tag "run_pepquery"

    container "proteomics/neoflow:latest"
    
    cpus "$params.cpu"

    publishDir "${out_dir}/pepquery/", mode: "copy", overwrite: true

    input:
    file(novel_peptide_tsv)
    file(pv_refdb)
    file(ms_data)

    output:
    file ("pepquery") into pepquery_res_folder

    script:
    if(ms_data.isFile()){
        // single file
        """
        #!/bin/sh
        cat 
        java -Xmx10g -jar /opt/pepquery-1.6.2/pepquery-1.6.2.jar \
              -ms ${ms_data} \
              -pep ${novel_peptide_tsv} \
              -db ${pv_refdb} \
              -fixMod ${pv_fixmod} \
              -varMod ${pv_varmod} \
              -cpu 0 \
              -minScore 12 \
              -tol ${pv_tol} \
              -tolu ${pv_tolu} \
              -itol ${pv_itol} \
              -n 10000 \
              -um \
              -m 1 \
              -prefix neoflow \
              -o pepquery
        """
    }else{
        // a folder
        """
        #!/bin/sh
        cat ${ms_data}/*.mgf > all_msms.mgf 
        java -Xmx10g -jar /opt/pepquery-1.6.2/pepquery-1.6.2.jar \
              -ms all_msms.mgf \
              -pep ${novel_peptide_tsv} \
              -db ${pv_refdb} \
              -fixMod ${pv_fixmod} \
              -varMod ${pv_varmod} \
              -cpu 0 \
              -minScore 12 \
              -tol ${pv_tol} \
              -tolu ${pv_tolu} \
              -itol ${pv_itol} \
              -n 10000 \
              -um \
              -m 1 \
              -prefix neoflow \
              -o pepquery
        """
        
    }
    
}

process add_pepquery_validation {

    tag "add_pepquery_validation"

    container "proteomics/pga:latest"
    
    cpus "$params.cpu"

    publishDir "${out_dir}/novel_peptide_identification/", mode: "copy", overwrite: true

    input:
    file novel_psm_tsv
    file pepquery_res_folder

    output:
    file("novel_peptides_psm_pepquery.tsv") into novel_peptides_psm_pepquery

    script:
    """
    #!/usr/bin/env /usr/local/bin/Rscript
    library(dplyr)
    library(readr)
    library(stringr)

    psms <- read_tsv("${novel_psm_tsv}")

    psm_rank_file = "${pepquery_res_folder}/psm_rank.txt"
    if(file.exists(psm_rank_file)){
        psm_rank <- read_tsv(psm_rank_file)
	if("n_ptm" %in% names(psm_rank)){
            psm_rank <- psm_rank %>% filter(pvalue<=0.01,n_ptm==0,rank==1)
            psms\$pepquery <- ifelse(psms\$peptide %in% psm_rank\$peptide,1,0)
	}else{
	    psms\$pepquery <- 0	
	}
    }else{
        psms\$pepquery <- 0
    }
    psms %>% write_tsv("novel_peptides_psm_pepquery.tsv")
    """

}


if (params.rt_validation) {
    
    process prepare_autort_data{

        tag "prepare_autort_data"

        container "proteomics/pga:latest"
	
	cpus "$params.cpu"

        input:
        file pga_re_dir from pga_result_folder_for_autort
        file(ms_data)
        file(pepquery_res_folder)
        output:
        file("autort_training_prediction_data") into autort_input_data

        script:
        """
        #!/bin/sh
        mkdir mgf_index
        java -cp ${baseDir}/bin/PDV-1.5.4_generate_index/PDV-1.5.4_generate_index.jar PDVGUI.generate_index ${ms_data} mgf_index/
        Rscript ${baseDir}/bin/rt_validation_data_preprocessing.R \
            psm_level/global_fdr/pga-peptideSummary.txt \
            peptide_level/global_fdr/pga-peptideSummary.txt \
            mgf_index/ \
            ${search_engine} \
            autort_training_prediction_data
        """
    }

    process run_autort{

        tag "run_autort"

        container "proteomics/autort:latest"
	
	cpus "$params.cpu"

        publishDir "${out_dir}/rt_validation/", mode: "copy", overwrite: true

        input:
        file(autort_input_data)

        output:
        //file(mgf_dir) into psm_fdr_ch
        file("rt_prediction") into final_res

        script:
        """
        #!/bin/sh
        python ${baseDir}/bin/run_autort.py ${autort_input_data} rt_prediction
        """
    }

    

}