nf_traelseq_rescue

#!/usr/bin/env nextflow

// Enable modules
nextflow.enable.dsl=2

params.outdir = "."
params.genome = ""
params.verbose = false
//params.single_end = false  // default mode is auto-detect. NOTE: params are handed over automatically 
params.single_end = true
params.fastqc_args = ''
params.trael_preprocessing_args = ''
params.fastq_screen_args = ''
params.trim_galore_args = ''
params.bowtie2_args = ''
params.samtools_sort_args = ''
params.samtools_index_args = ''
params.umibam2_args = ''
params.index = '123456789' // we want to supply the indexes to use
params.no_index = false
params.multiqc_args = ''
params.no_output = false

params.help = false
// Show help message and exit
if (params.help){
    helpMessage()
    exit 0
}

params.list_genomes = false;
if (params.list_genomes){
    println ("[WORKFLOW] List genomes selected")
}

if (params.verbose){
    println ("[WORKFLOW] FASTQC ARGS: "           + params.fastqc_args)
    println ("[WORKFLOW] TrAEL-PREPROCESSING ARGS: " + params.trael_preprocessing_args)
    println ("[WORKFLOW] FASTQ SCREEN ARGS ARE: " + params.fastq_screen_args)
    println ("[WORKFLOW] TRIM GALORE ARGS: "      + params.trim_galore_args)
    println ("[WORKFLOW] BOWTIE2 ARGS: "          + params.bowtie2_args)
    println ("[WORKFLOW] UMIBAM ARGS ARE: "       + params.umibam2_args)
    println ("[WORKFLOW] MULTIQC ARGS: "          + params.multiqc_args)
}


include { makeFilesChannel; getFileBaseNames } from './nf_modules/files.mod.nf'
include { getGenome }                          from './nf_modules/genomes.mod.nf'
include { listGenomes }                        from './nf_modules/genomes.mod.nf'

if (params.list_genomes){
    listGenomes()  // this lists all available genomes, and exits
}
genome = getGenome(params.genome)

include { TRAEL_PREPROCESSING_INDEXING } from './nf_modules/trael_preprocessing.mod.nf'
include { TRAEL_PREPROCESSING } from './nf_modules/trael_preprocessing.mod.nf'
include { SORT_TRAEL_NAMES }    from './nf_modules/trael_preprocessing.mod.nf'
include { FASTQC }              from './nf_modules/fastqc.mod.nf'    params(nogroup: true)
include { FASTQC as FASTQC2 }   from './nf_modules/fastqc.mod.nf'    params(nogroup: true)
include { FASTQC as FASTQC3 }   from './nf_modules/fastqc.mod.nf'    params(nogroup: true)
include { FASTQ_SCREEN }        from './nf_modules/fastq_screen.mod.nf'
include { TRIM_GALORE }         from './nf_modules/trim_galore.mod.nf' params(single_end: '--single_end')
include { SAMTOOLS_SORT }       from './nf_modules/samtools.mod.nf'  
include { SAMTOOLS_INDEX2 }      from './nf_modules/samtools.mod.nf' 
include { BOWTIE2 }             from './nf_modules/bowtie2.mod.nf'   params(genome: genome, no_output: params.no_output)
// params(genome: genome, local: "--local")
include { UMIBAM2 }              from './nf_modules/umibam2.mod.nf' // new umibam needs putting into nextflow pipeline
include { MULTIQC }              from './nf_modules/multiqc.mod.nf' 
// https://bioinformaticsworkbook.org/dataAnalysis/nextflow/02_creatingAworkflow.html#gsc.tab=0 

file_ch = makeFilesChannel(args)

matcher_string = ".*index[" + params.index + "].*"

/* 
At first this seemed like complete overkill but I couldn't work out a nextflow way of doing this
The fq.gz files are "not empty" so can't use ifEmpty. Bash also doesn't recognise them 
as empty so this seems like a reasonable way to check. 
The line below can be used to filter the output.
filt = CHECK_FQ_GZ_NOT_EMPTY.out.reads.filter { it.first().countLines() > 0 }
*/
process CHECK_FQ_GZ_NOT_EMPTY {

    input:
    path(reads)

    output:
    tuple path ("*fq"), path (reads), emit: reads

    script:
    """
    zcat ${reads} | head -n 5 > first_line.fq
    """
}

// same approach as the CHECK_FQ_GZ_NOT_EMPTY process above but for bam files
process CHECK_BAMFILE_NOT_EMPTY {

    input:
    path(bam)

    output:
    tuple path ("*sam"), path (bam), emit: bam

    script:
    """
    module load samtools
    samtools view ${bam} | head -n 1 > first_line.sam
    """
}

process TRAEL_PREPROCESSING_INDEXING_RESCUE {
	
	tag "$name" // Adds name to job submission instead of (1), (2) etc.
		
    input:
	    tuple val(name), path(reads)
		val (outputdir)
		val (trael_preprocessing_args)
		val (verbose)

	output:
		path "*.txt", optional: true, emit: stats 
        path ("*UMIed*.fastq.gz"), emit: reads

	publishDir "$outputdir",
		mode: "link", overwrite: true, enabled: !params.no_output

	script:
		if (verbose){
			println ("[MODULE] TrAEL-PREPROCESSING ARGS: " + trael_preprocessing_args)
		}

		/* Run the TrAELseq_preprocessing_UMIplusBarcode script
			This splits the input fastq file by the 9 TrAEL barcodes and T or no T at position 13, 
			so produces 20 output fastq files (9 + unassigned)*2 for T and noT. 
		*/
		"""
		module load python
		/bi/group/bioinf/Jon_Houseley/trael_rescue_lane8526/TrAELseq_UMIplusBarcode_rescue_lane8526.py	${reads}
        #/bi/group/bioinf/Jon_Houseley/trael_rerescue/sample5693_lane8149_8160/TrAELseq_preprocessing_UMIplusBarcode_rerescue_lane8149_8160.py	${reads}
		"""

}


workflow {

    main:
        //FASTQC          (file_ch, params.outdir, params.fastqc_args, params.verbose)
        //FASTQ_SCREEN    (file_ch, params.outdir, params.fastq_screen_args, params.verbose)

        if ( params.no_index ) {
            TRAEL_PREPROCESSING  (file_ch, params.outdir, params.trael_preprocessing_args, params.verbose)
            CHECK_FQ_GZ_NOT_EMPTY(TRAEL_PREPROCESSING.out.reads.flatten())

        }
        else  {
            TRAEL_PREPROCESSING_INDEXING_RESCUE  (file_ch, params.outdir, params.trael_preprocessing_args, params.verbose)
            CHECK_FQ_GZ_NOT_EMPTY(TRAEL_PREPROCESSING_INDEXING_RESCUE.out.reads.flatten())
        }       

        /* The TRAEL_PREPROCESSING scripts produce a set of fastq files for each name, but some may be empty. 
        We need to filter these out as some processes e.g. trim_galore will fail if empty files are passed. 
        Also, we don't want to keep processing a load of empty files. 
        */
                      
        filt = CHECK_FQ_GZ_NOT_EMPTY.out.reads.filter { it.first().countLines() > 0 }.map { it.last() }

        // or in 2 steps
        //filt = CHECK_FQ_GZ_NOT_EMPTY.out.reads.filter { it.first().countLines() > 0 }
        //SORT_TRAEL_NAMES(filt.map { it.last() })

        // Required input for bowtie2 module is a tuple of one name, one path, not just a set of paths
        // which is what TRAEL_PREPROCESSING_INDEXING produces
        SORT_TRAEL_NAMES(filt)


        if ( params.no_index ) {
            trael_ch = SORT_TRAEL_NAMES.out.reads
        }
        else {
            // only keep the selected indexes
            trael_ch = SORT_TRAEL_NAMES.out.reads.filter( ~ /${matcher_string}/ )
        }
        //FASTQC2     (trael_ch, params.outdir, params.fastqc_args, params.verbose)
        TRIM_GALORE (trael_ch, params.outdir, params.trim_galore_args, params.verbose)
        FASTQC3     (TRIM_GALORE.out.reads, params.outdir, params.fastqc_args, params.verbose)
        BOWTIE2     (TRIM_GALORE.out.reads, params.outdir, params.bowtie2_args, params.verbose)

        // we just need the bam files for SAMTOOLS_SORT and INDEX, not a tuple of the name and the bam
        mapped_bams = BOWTIE2.out.bam.map { it.last() }

        // need to filter out any empty bam files from the channel - umibam fails if passed an empty bam file
        CHECK_BAMFILE_NOT_EMPTY(mapped_bams)
        filt_sam_bam = CHECK_BAMFILE_NOT_EMPTY.out.bam.filter { it.first().countLines() > 0 }
        filt_bam = filt_sam_bam.map { it.last() }
        SAMTOOLS_SORT  (filt_bam, params.outdir, params.samtools_sort_args, params.verbose)
        //SAMTOOLS_INDEX (SAMTOOLS_SORT.out.bam, params.outdir, params.samtools_index_args, params.verbose)
      				        
        // Deduplication
        //UMIBAM2 (SAMTOOLS_SORT.out.bam, SAMTOOLS_INDEX.out.bai, params.outdir, params.umibam2_args, params.verbose)

        SAMTOOLS_INDEX2 (SAMTOOLS_SORT.out.bam, params.outdir, params.samtools_index_args, params.verbose)
      				        
        // Deduplication
        //UMIBAM2 (SAMTOOLS_SORT.out.bam, SAMTOOLS_INDEX.out.bai, params.outdir, params.umibam2_args, params.verbose)
        UMIBAM2 (SAMTOOLS_INDEX2.out.bam, SAMTOOLS_INDEX2.out.bai, params.outdir, params.umibam2_args, params.verbose)


        // merging channels for MultiQC
      /*  multiqc_ch = FASTQC.out.report.mix(
            FASTQ_SCREEN.out.report.ifEmpty([]),
            FASTQC2.out.report.ifEmpty([]),
            FASTQC3.out.report.ifEmpty([]),
            TRIM_GALORE.out.report.ifEmpty([]),
            BOWTIE2.out.stats.ifEmpty([]),
            UMIBAM2.out.report.ifEmpty([]),
        ).collect()
        */
        //multiqc_ch.subscribe {  println "Got: $it"  }

        multiqc_ch = FASTQC3.out.report.mix( 
            TRIM_GALORE.out.report.ifEmpty([]),
            BOWTIE2.out.stats.ifEmpty([]),
            UMIBAM2.out.report.ifEmpty([]),
        ).collect()

        MULTIQC  (multiqc_ch, params.outdir, params.multiqc_args, params.verbose)
        
}

// Since workflows with very long command lines tend to fail to get rendered at all, I was experimenting with a
// minimal execution summary report so we at least know what the working directory was...
workflow.onComplete {

    def msg = """\
        Pipeline execution summary
        ---------------------------
        Jobname     : ${workflow.runName}
        Completed at: ${workflow.complete}
        Duration    : ${workflow.duration}
        Success     : ${workflow.success}
        workDir     : ${workflow.workDir}
        exit status : ${workflow.exitStatus}
        """
    .stripIndent()

    sendMail(to: "${workflow.userName}@babraham.ac.uk", subject: 'Minimal pipeline execution report', body: msg)
}

def helpMessage() {
 
    log.info"""
    >>

    SYNOPSIS:

    This workflow runs a TrAEL-seq processing pipeline from FastQ files, including QC, contamination QC,
    splitting by TrAEL barcodes, quality-/adapter trimming, alignments to a genome using Bowtie 2, 
    deduplication based on position and UMI, and finally generates an aggregate QC report. 

    Here is a graphical representation of the workflow:

    --- FastQC
    --- FastQ Screen
    --- TrAEL preprocessing 
        | 
        --- FastQC
        --- Trim Galore
            |
            --- FastQC
            --- Bowtie2  
                |
                --- Samtools sort 
                --- Samtools index  
                    |
                    umibam2
    --- MultiQC*
        
    * This step runs only once ALL other jobs have completed.
     
    By default all these steps are submitted as jobs to the Babraham stone compute cluster.

    By default, the involved tools are run in the following way:
    ------------------------------------------------------------
    FastQC:         defaults (-q)
    FastQ Screen:   defaults (Bowtie 2; local mode)
    Trim Galore:    defaults (adapter auto-detection)
    Bowtie 2:       end-to-end mode; '--no-unal'; for paired-end files: '--no-mixed --no-discordant' (concordant PE alignmnents only)
                    
    To add additional parameters to any of the programs, consider supplying tool-specific arguments (see --toolname_args="..." below).

            ==============================================================================================================

            
    USAGE:
    
    nf_traelseq [options] --genome <genomeID> --index <TrAEL indexes> <input files>
    
    Mandatory arguments:
    ====================

      <input files>                   List of input files, e.g. '*fastq.gz' or '*fq.gz'. Files are automatically processed as
                                      single-end (SE) or paired end (PE) files (if file pairs share the same base-name, and differ only
                                      by a read number, e.g. 'base_name_R1.fastq.gz' and 'base_name_R2.fastq.gz' (or R3, R4). For
                                      PE files, only Read 1 is run through FastQ Screen (as typically R1 and R2 produce nearly identical
                                      contamination profiles). To run PE files in single-end mode, please see '--single_end' below.

      --genome [str]                  Genome build ID to be used for the alignment, e.g. GRCh38 (latest human genome) or GRCm38
                                      (latest mouse genome build). To list all available genomes, see '--list_genomes' below.


    Optional arguments:
    ====================

      --index [int]                  TrAEL indexes that are expected to be present. Provide the index numbers without spaces, 
                                        e.g. --index 245 for indexes 2,4 and 5.  [Default: all indexes, equivalent to --index 123456789]

                                     These are the TrAEL indexes:   
                                        barcode_1 = "AGTC" 
                                        barcode_2 = "GACT" 
                                        barcode_3 = "CTTG"
                                        barcode_4 = "TCGA"
                                        barcode_5 = "AAGG"
                                        barcode_6 = "TTCC"
                                        barcode_7 = "GTGC"
                                        barcode_8 = "GCCA"
                                        barcode_9 = "GATG"
                                  

    Tool-specific options:
    ======================

    For all following options, please note that the format: ="your options" needs to be strictly adhered to in order to work correctly.

      --fastqc_args="[str]"           This option can take any number of options that are compatible with FastQC to modify its default
                                      behaviour. For more detailed information on available options please refer to the FastQC documentation,
                                      or run 'fastqc --help' on the command line. As an example, to run FastQC without grouping of bases when
                                      reads are >50bp and use a specific file with non-default adapter sequences, use: 
                                      ' --fastqc_args="--nogroup --adapters ./non_default_adapter_file.txt" '. [Default: None]

      --fastq_screen_args="[str]"     This option can take any number of options that are compatible with FastQ Screen to modify its
                                      default behaviour. For more detailed information on available options please refer to the FastQ Screen
                                      documentation, or run 'fastq_screen --help' on the command line. For instance, to process a bisulfite
                                      converted library with fairly relaxed parameters, you could use: 
                                      ' --fastq_screen_args="--bisulfite --score_min L,0,-0.6" '. [Default: None]

      --trim_galore_args="[str]"      This option can take any number of options that are compatible with Trim Galore to modify its
                                      default trimming behaviour. For more detailed information on available options please refer
                                      to the Trim Galore User Guide, or run 'trim_galore --help' on the command line. As an example, to trim
                                      off the first 10bp from the 5' of R1 and 5bp of R2, use: 
                                      ' --trim_galore_args="--clip_r1 10 --clip_r2 5" '. [Default: None]

    
      --bowtie2_args="[str]"          This option can take any number of options that are compatible with Bowtie 2 to modify its
                                      default mapping behaviour. For more detailed information on available options please refer
                                      to the Bowtie 2 User Guide, or run 'bowtie2 --help' on the command line. As an example, to
                                      run somewhat more stringent alignments for only 1 million sequences, use:
                                      ' --bowtie2_args="-u 1000000 --score_min L,0,-0.4" '. [Default: None]

    Other options:
    ==============

      --outdir [str]                  Path to the output directory. [Default: current working directory]

      --list_genomes                  List all genome builds that are currently available to choose from. To see this list
                                      of available genomes with more detailed information about paths and indexes, run
                                      the command as '--list_genomes --verbose'
    
      --single_end                    Force files of a read pair to be treated as single-end files. [Default: auto-detect]
      
      --verbose                       More verbose status messages. [Default: OFF]
      --help                          Displays this help message and exits.

    Workflow options:
    =================

    Please note the single '-' hyphen for the following options!

      -resume                         If a pipeline workflow has been interrupted or stopped (e.g. by accidentally closing a laptop),
                                      this option will attempt to resume the workflow at the point it got interrupted by using
                                      Nextflow's caching mechanism. This may save a lot of time.

      -bg                             Sends the entire workflow into the background, thus disconnecting it from the terminal session.
                                      This option launches a daemon process (which will keep running on the headnode) that watches over
                                      your workflow, and submits new jobs to the SLURM queue as required. Use this option for big pipeline
                                      jobs, or whenever you do not want to watch the status progress yourself. Upon completion, the
                                      pipeline will send you an email with the job details. This option is HIGHLY RECOMMENDED!

      -process.executor=local         Temporarily changes where the workflow is executed to the 'local' machine. See also Nextflow config
                                      file for more details. [Default: slurm] 
    

    <<
    """.stripIndent()

}