# Example code for FAQ Good practices: “What is the best way of handling data in file pairs (or triplets etc.)?”

> See also

> - @collate

>     #!/usr/bin/env python
>     import sys, os

>     from ruffus import *
>     import ruffus.cmdline as cmdline
>     from subprocess import check_call

>     parser = cmdline.get_argparse(description="Parimala's pipeline?")

>     #                                                                                 .
>     #   Very flexible handling of input files                                         .
>     #                                                                                 .
>     #      input files can be specified flexibly as:                                  .
>     #                 --input a.fastq b.fastq                                         .
>     #                 --input a.fastq --input b.fastq                                 .
>     #                 --input *.fastq --input other/*.fastq                           .
>     #                 --input "*.fastq"                                               .
>     #                                                                                 .
>     #       The last form is expanded in the script and avoids limitations on command .
>     #           line lengths                                                          .
>     #                                                                                 .
>     parser.add_argument('-i', '--input', nargs='+', metavar="FILE", action="append", help = "Fastq files")

>     options = parser.parse_args()

>     #  standard python logger which can be synchronised across concurrent Ruffus tasks
>     logger, logger_mutex = cmdline.setup_logging ("PARIMALA", options.log_file, options.verbose)

>     #                                                                                .
>     #   Useful code to turn input files into a flat list                             .
>     #                                                                                .
>     from glob import glob
>     original_data_files = [fn for grouped in options.input for glob_spec in grouped for fn in glob(glob_spec)] if options.input else []
>     if not original_data_files:
>         original_data_files = [["C1W1_R1.fastq.gz", "C1W1_R2.fastq.gz"]]
>         #raise Exception ("No matching files specified with --input.")

>     #   <<<----  pipelined functions go here

>     #_________________________________________________________________________________
>     #                                                                                .
>     #   Group together file pairs                                                    .
>     #_________________________________________________________________________________
>     @collate(original_data_files,
>                 # match file name up to the "R1.fastq.gz"
>                 formatter("([^/]+)R[12].fastq.gz$"),
>                 # Create output parameter supplied to next task
>                 ["{path[0]}/{1[0]}paired.R1.fastq.gz",  # paired file 1
>                  "{path[0]}/{1[0]}paired.R2.fastq.gz"], # paired file 2
>                 # Extra parameters for our own convenience and use
>                 ["{path[0]}/{1[0]}unpaired.R1.fastq.gz",  # unpaired file 1
>                  "{path[0]}/{1[0]}unpaired.R2.fastq.gz"], # unpaired file 2
>                 logger, logger_mutex)
>     def trim_fastq(input_files, output_paired_files, discarded_unpaired_files, logger, logger_mutex):
>         if len(input_files) != 2:
>             raise Exception("One of read pairs %s missing" % (input_files,))
>         cmd = ("java -jar ~/SPRING-SUMMER_2014/Softwares/Trimmomatic/Trimmomatic-0.32/trimmomatic-0.32.jar "
>                " PE -phred33  "
>                " {input_files[0]} {input_files[1]} "
>                " {output_paired_files[0]} {output_paired_files[1]} "
>                " {discarded_unpaired_files[0]} {discarded_unpaired_files[1]} "
>                " LEADING:30 TRAILING:30 SLIDINGWINDOW:4:15 MINLEN:50 "
>             )

>         check_call(cmd.format(**locals()))

>         with logger_mutex:
>                 logger.debug("Hooray trim_fastq worked")

>     #_________________________________________________________________________________
>     #                                                                                .
>     #   Each file pair now makes its way down the rest of the pipeline as            .
>     #       a couple                                                                 .
>     #_________________________________________________________________________________
>     @transform(trim_fastq,
>                 # regular expression match on first of pe files
>                 formatter("([^/]+)paired.R1.fastq.gz$"),
>                 # Output parameter supplied to next task
>                 "{path[0]}/{1[0]}.sam"

>                # Extra parameters for our own convenience and use
>                 "{path[0]}/{1[0]}.pe_soap_pe",        # soap intermediate file
>                 "{path[0]}/{1[0]}.pe_soap_se",        # soap intermediate file
>                 logger, logger_mutex)
>     def align_seq(input_files, output_file, soap_pe_output_file, soap_se_output_file, logger, logger_mutex):
>         if len(input_files) != 2:
>             raise Exception("One of read pairs %s missing" % (input_files,))
>         cmd = ("~/SPRING-SUMMER_2014/Softwares/soap2.21release/soap "
>                 " -a {input_files[0]} "
>                 " -b {input_files[1]} "
>                 " -D Y55_genome.fa.index* "
>                 " -o {soap_pe_output_file} -2 {soap_se_output_file} -m 400 -x 600")

>         check_call(cmd.format(**locals()))

>         #Soap_to_sam
>         cmd = " perl ~/SPRING-SUMMER_2014/Softwares/soap2sam.pl -p {soap_pe_output_file} > {output_file}"

>         check_call(cmd.format(**locals()))

>         with logger_mutex:
>                 logger.debug("Hooray align_seq worked")

>     cmdline.run (options)