##############################################################################
#
#   Kennedy Institute of Rheumatology
#
#   $Id$
#
#   Copyright (C) 2015 <NAME>
#
#   This program is free software; you can redistribute it and/or
#   modify it under the terms of the GNU General Public License
#   as published by the Free Software Foundation; either version 2
#   of the License, or (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
###############################################################################

"""
===========================
Pipeline cram2fastq
===========================

:Author: <NAME>
:Release: $Id$
:Date: |today|
:Tags: Python

Overview
========

This pipeline coverts Sanger CRAM files to fastq.gz,
optionally quality trimming and reconciling the fastq files

Usage
=====

See :ref:`PipelineSettingUp` and :ref:`PipelineRunning` on general
information how to use CGAT pipelines.

Configuration
-------------

The pipeline requires a configured :file:`pipeline.ini` file.
CGATReport report requires a :file:`conf.py` and optionally a
:file:`cgatreport.ini` file (see :ref:`PipelineReporting`).

Default configuration files can be generated by executing:

   python <srcdir>/pipeline_cram2fastq.py config

Input files
-----------

Requirements
------------


On top of the default CGAT setup, the pipeline requires the following
software to be in the path:

.. Add any additional external requirements such as 3rd party software
   or R modules below:

Requirements:

* samtools >= 1.1

Pipeline output
===============


Glossary
========


Code
====

"""

from ruffus import *

import sys
import os
import glob
import sqlite3

from CGATCore import Experiment as E
from CGATCore import Pipeline as P
from CGATCore import Database as DB

import pysam

# -------------------------- < parse parameters > --------------------------- #

# load options from the config file
PARAMS = P.get_parameters(
    ["%s/pipeline.yml" % os.path.splitext(__file__)[0],
     "../pipeline.yml",
     "pipeline.yml"])


# ----------------------- < pipeline configuration > ------------------------ #

if len(sys.argv) > 1:
    if(sys.argv[1] == "config") and __name__ == "__main__":
        sys.exit(P.main(sys.argv))


# ------------------------< specific pipeline tasks >------------------------ #

@follows(mkdir("validate.cram.dir"))
@transform(glob.glob("data.dir/*.cram"),
           regex(r".*/(.*).cram"),
           [r"validate.cram.dir/\1.validate", r"validate.cram.dir/\1.quality"])
def validateCramFiles(infile, outfiles):
    '''Validate CRAM files by exit status of
       cramtools qstat. Save the quality scores of cram files.
    '''

    outfile, outfile_quality = outfiles

    statement = '''temp_quality=`mktemp -p %(cluster_tmpdir)s`;
                   cramtools qstat -I %(infile)s > $temp_quality;
                   echo $? > %(outfile)s;
                   cat $temp_quality
                   | awk '{OFS="\\t"} {print $1,$2}'
                   > %(outfile_quality)s;
                   rm $temp_quality;
                '''

    P.run(statement)


@follows(validateCramFiles)
@merge(validateCramFiles,
       "validate.cram.dir/summary.txt")
def inspectValidations(infiles, outfile):
    '''Check that all crams pass validation or
       raise an Error.'''

    validation_files = [fn
                        for filenames in infiles
                        for fn in filenames
                        if fn.endswith(".validate")]

    outfile_handle = open(outfile, "w")

    exit_states = []
    for validation_file in validation_files:

        with open(validation_file, "r") as vf_handle:
            exit_status = vf_handle.read().strip("\n")

        exit_states.append(int(exit_status))
        outfile_handle.write("\t".join([validation_file, exit_status])+"\n")

    outfile_handle.close()

    if sum(exit_states) != 0:
        raise ValueError("One or more cram files failed validation")


@follows(validateCramFiles)
@merge(validateCramFiles,
       "validate.cram.dir/cram_quality.load")
def loadCramQuality(infiles, outfile):
    ''' Load the quality scores for the different cells
        into the database (summarized table).
    '''

    quality_files = [fn
                     for filenames in infiles
                     for fn in filenames
                     if fn.endswith(".quality")]

    P.concatenate_and_load(quality_files, outfile,
                           regex_filename="validate.cram.dir/(.*).quality",
                           cat="track",
                           has_titles=False,
                           header="cramID,number_reads,cram_quality_score")


@follows(inspectValidations,
         mkdir("cell.info.dir"))
@merge(glob.glob("data.dir/*.cram"),
       "cell.info.dir/cells.txt")
def extractSampleInformation(infiles, outfile):
    '''Make a table of cells and corresponding cram files'''

    # build a dictionary of cell to cram file mappings
    cells = {}
    for cram_file in infiles:
        cram = pysam.AlignmentFile(cram_file, "rc")
        print(cram.header)
        cell = cram.header["RG"][0]["SM"]
        if cell not in cells.keys():
            cells[cell] = [cram_file]
        else:
            cells[cell].append(cram_file)
        cram.close()

    # write out a per-cell list of cram files
    outdir = os.path.dirname(outfile)

    outfile_handle = open(outfile, "w")
    outfile_handle.write("#cell\tcram_files\n")

    for cell in cells.keys():
        outfile_handle.write("\t".join([cell, ",".join(cells[cell])])+"\n")

    outfile_handle.close()


@split(extractSampleInformation,
       "cell.info.dir/*.cell")
def cellCramLists(infile, outfiles):
    '''Make a per-cell file containing the cram file(s)
       corresponding to the cell'''

    out_dir = os.path.dirname(infile)

    with open(infile, "r") as cell_list:
        for record in cell_list:
            if record.startswith("#"):
                continue
            cell, cram_list = record.strip("\n").split("\t")
            crams = cram_list.split(",")

            cell_outfile_name = os.path.join(out_dir, cell+".cell")
            with open(cell_outfile_name, "w") as cell_file_handle:
                for cram in crams:
                    cell_file_handle.write(cram+"\n")


@follows(mkdir("fastq.dir"),
         mkdir("fastq.temp.dir"),
         extractSampleInformation)
@transform(cellCramLists,
           regex(r".*/(.*).cell"),
           (r"fastq.dir/\1.fastq.1.gz",
            r"fastq.dir/\1.fastq.2.gz"))
def cram2fastq(infile, outfiles):
    '''Convert Sanger CRAM files to Fastq format
       Takes care of merging, quality trimming
       and pair reconciliation.
       Intermediate files are not kept by default.'''

    # TODO: make quality trimming optional.

    ###################################
    # set variables and open a log file
    ###################################

    cell_name = os.path.basename(infile)[:-len(".cell")]
    out_dir = os.path.dirname(outfiles[0])
    temp_dir = "fastq.temp.dir"

    log_file = os.path.join(temp_dir,
                            cell_name + ".fastq.extraction.log")

    log = open(log_file, "w")
    log.write("Fastq extraction log file for %(infile)s\n\n")

    def _merge_dicts(a, b):
        x = a.copy()
        x.update(b)
        return(x)

    temp_files = []

    # ##############################################
    # Extract per-end Fastq(s) from the cram file(s)
    # ##############################################

    raw_fastq_names = []
    with open(infile, "rb") as cram_files:

        for line in cram_files:

            cram = line.strip()
            cram_basename = os.path.basename(cram)[:-len(".cram")]
            raw_fastq_name = os.path.join(temp_dir, cram_basename)
            raw_fastq_names.append(raw_fastq_name)

            job_memory = PARAMS["preprocess_memory"]
            statement = '''cramtools fastq
                                       --enumerate
                                       --reverse
                                       -F %(raw_fastq_name)s
                                       -I %(cram)s
                                       --gzip
                         '''
            log.write("Extracting fastqs from %(cram)s:" % locals() + "\n")
            log.write(statement % locals() + "\n")

            P.run(statement)

            log.write("done.\n\n")

    # ####################################
    # Perform quality trimming
    # Merging is also taken care of here.
    # ####################################

    quality = PARAMS["preprocess_quality_threshold"]
    minlen = PARAMS["preprocess_min_length"]
    trim = PARAMS["preprocess_trim"]

    trimmed_fastq_prefix = os.path.join(temp_dir, cell_name)

    trimmed_fastq_files = []

    # fastq(s) for each end are quality trimmed separately
    for end in ["_1", "_2"]:

        raw_fastqs = [x + end + ".fastq.gz" for x in raw_fastq_names]
        temp_files += raw_fastqs

        fastq_list = " ".join(raw_fastqs)

        trimmed_fastq_name = trimmed_fastq_prefix + end + ".trimmed.fastq.gz"
        trimmed_fastq_files.append(trimmed_fastq_name)

        log.write(">> Quality trimming %(fastq_list)s: " % locals() + "\n")
        if trim:
            statement = '''zcat %(fastq_list)s
                       | fastq_quality_trimmer
                           -Q33
                           -t %(quality)s
                           -l %(minlen)s
                       | gzip -c
                       > %(trimmed_fastq_name)s
                       '''
        else:
            statement = '''zcat %(fastq_list)s
                       | gzip -c
                       > %(trimmed_fastq_name)s
                       '''
        log.write(statement % _merge_dicts(PARAMS, locals()) + "\n")
        P.run(statement)
        log.write("done. \n\n")

    # ##################
    # Reconcile the ends
    # ##################

    if PARAMS["preprocess_reconcile"] != "False":

        temp_files += trimmed_fastq_files

        end1, end2 = trimmed_fastq_files

        reconciled_fastq_prefix = outfiles[0][:-len(".1.gz")]

        log.write(">> Reconciling pairs, %(end1)s & %(end2)s: "
                  % locals() + "\n")
        statement = '''python %(scriptsdir)s/fastqs2fastqs.py
                       %(end1)s %(end2)s
                       --method reconcile
                       --chop
                       --unpaired
                       -o "%(reconciled_fastq_prefix)s.%%s.gz";
                    '''
        log.write(statement % _merge_dicts(PARAMS, locals()) + "\n")
        P.run(statement)
        log.write("done\n\n")

    else:
        trimmed_fastq_prefix = outfiles[0][:-len(".1.gz")]
        for end in trimmed_fastq_files:
            if "1.trimmed" in end:
                endn = "1"
            else:
                endn = "2"

            trimmed_end_name = ".".join([trimmed_fastq_prefix, endn, "gz"])
            os.symlink(os.path.abspath(end), trimmed_end_name)

    ##############################
    # Clean up the temporary files
    ##############################

    if PARAMS["keep_temporary"] == 0:

        temp_file_list = " ".join(temp_files)

        # record files sizes and md5 checksums of the temporary files
        log.write(">> Recording sizes and checksums of temporary files:\n")
        statement = '''ls -l %(temp_file_list)s
                       > %(temp_dir)s/%(cell_name)s.ls;
                       checkpoint;
                       md5sum %(temp_file_list)s
                       > %(temp_dir)s/%(cell_name)s.md5;
                    '''
        log.write(statement % locals() + "\n")
        P.run(statement)
        log.write("done\n\n")

        # unlink (delete) the temporary files
        log.write(">> unlinking temporary files: " + temp_file_list + "\n")

        for temp_file in temp_files:
            os.unlink(temp_file)

        log.write("tempororay files unlinked\n")

    log.close()


# ---------------------< generic pipeline tasks >---------------------------- #

@follows(cram2fastq, loadCramQuality)
def full():
    pass


# ########################################################################### #

if __name__ == "__main__":
    sys.exit(P.main(sys.argv))

##############################################################################
#
#   Kennedy Institute of Rheumatology
#
#   $Id$
#
#   Copyright (C) 2015 <NAME>
#
#   This program is free software; you can redistribute it and/or
#   modify it under the terms of the GNU General Public License
#   as published by the Free Software Foundation; either version 2
#   of the License, or (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
###############################################################################

"""
===========================
Pipeline cram2fastq
===========================

:Author: <NAME>
:Release: $Id$
:Date: |today|
:Tags: Python

Overview
========

This pipeline coverts Sanger CRAM files to fastq.gz,
optionally quality trimming and reconciling the fastq files

Usage
=====

See :ref:`PipelineSettingUp` and :ref:`PipelineRunning` on general
information how to use CGAT pipelines.

Configuration
-------------

The pipeline requires a configured :file:`pipeline.ini` file.
CGATReport report requires a :file:`conf.py` and optionally a
:file:`cgatreport.ini` file (see :ref:`PipelineReporting`).

Default configuration files can be generated by executing:

   python <srcdir>/pipeline_cram2fastq.py config

Input files
-----------

Requirements
------------


On top of the default CGAT setup, the pipeline requires the following
software to be in the path:

.. Add any additional external requirements such as 3rd party software
   or R modules below:

Requirements:

* samtools >= 1.1

Pipeline output
===============


Glossary
========


Code
====

"""

from ruffus import *

import sys
import os
import glob
import sqlite3

from CGATCore import Experiment as E
from CGATCore import Pipeline as P
from CGATCore import Database as DB

import pysam

# -------------------------- < parse parameters > --------------------------- #

# load options from the config file
PARAMS = P.get_parameters(
    ["%s/pipeline.yml" % os.path.splitext(__file__)[0],
     "../pipeline.yml",
     "pipeline.yml"])


# ----------------------- < pipeline configuration > ------------------------ #

if len(sys.argv) > 1:
    if(sys.argv[1] == "config") and __name__ == "__main__":
        sys.exit(P.main(sys.argv))


# ------------------------< specific pipeline tasks >------------------------ #

@follows(mkdir("validate.cram.dir"))
@transform(glob.glob("data.dir/*.cram"),
           regex(r".*/(.*).cram"),
           [r"validate.cram.dir/\1.validate", r"validate.cram.dir/\1.quality"])
def validateCramFiles(infile, outfiles):
    '''Validate CRAM files by exit status of
       cramtools qstat. Save the quality scores of cram files.
    '''

    outfile, outfile_quality = outfiles

    statement = '''temp_quality=`mktemp -p %(cluster_tmpdir)s`;
                   cramtools qstat -I %(infile)s > $temp_quality;
                   echo $? > %(outfile)s;
                   cat $temp_quality
                   | awk '{OFS="\\t"} {print $1,$2}'
                   > %(outfile_quality)s;
                   rm $temp_quality;
                '''

    P.run(statement)


@follows(validateCramFiles)
@merge(validateCramFiles,
       "validate.cram.dir/summary.txt")
def inspectValidations(infiles, outfile):
    '''Check that all crams pass validation or
       raise an Error.'''

    validation_files = [fn
                        for filenames in infiles
                        for fn in filenames
                        if fn.endswith(".validate")]

    outfile_handle = open(outfile, "w")

    exit_states = []
    for validation_file in validation_files:

        with open(validation_file, "r") as vf_handle:
            exit_status = vf_handle.read().strip("\n")

        exit_states.append(int(exit_status))
        outfile_handle.write("\t".join([validation_file, exit_status])+"\n")

    outfile_handle.close()

    if sum(exit_states) != 0:
        raise ValueError("One or more cram files failed validation")


@follows(validateCramFiles)
@merge(validateCramFiles,
       "validate.cram.dir/cram_quality.load")
def loadCramQuality(infiles, outfile):
    ''' Load the quality scores for the different cells
        into the database (summarized table).
    '''

    quality_files = [fn
                     for filenames in infiles
                     for fn in filenames
                     if fn.endswith(".quality")]

    P.concatenate_and_load(quality_files, outfile,
                           regex_filename="validate.cram.dir/(.*).quality",
                           cat="track",
                           has_titles=False,
                           header="cramID,number_reads,cram_quality_score")


@follows(inspectValidations,
         mkdir("cell.info.dir"))
@merge(glob.glob("data.dir/*.cram"),
       "cell.info.dir/cells.txt")
def extractSampleInformation(infiles, outfile):
    '''Make a table of cells and corresponding cram files'''

    # build a dictionary of cell to cram file mappings
    cells = {}
    for cram_file in infiles:
        cram = pysam.AlignmentFile(cram_file, "rc")
        print(cram.header)
        cell = cram.header["RG"][0]["SM"]
        if cell not in cells.keys():
            cells[cell] = [cram_file]
        else:
            cells[cell].append(cram_file)
        cram.close()

    # write out a per-cell list of cram files
    outdir = os.path.dirname(outfile)

    outfile_handle = open(outfile, "w")
    outfile_handle.write("#cell\tcram_files\n")

    for cell in cells.keys():
        outfile_handle.write("\t".join([cell, ",".join(cells[cell])])+"\n")

    outfile_handle.close()


@split(extractSampleInformation,
       "cell.info.dir/*.cell")
def cellCramLists(infile, outfiles):
    '''Make a per-cell file containing the cram file(s)
       corresponding to the cell'''

    out_dir = os.path.dirname(infile)

    with open(infile, "r") as cell_list:
        for record in cell_list:
            if record.startswith("#"):
                continue
            cell, cram_list = record.strip("\n").split("\t")
            crams = cram_list.split(",")

            cell_outfile_name = os.path.join(out_dir, cell+".cell")
            with open(cell_outfile_name, "w") as cell_file_handle:
                for cram in crams:
                    cell_file_handle.write(cram+"\n")


@follows(mkdir("fastq.dir"),
         mkdir("fastq.temp.dir"),
         extractSampleInformation)
@transform(cellCramLists,
           regex(r".*/(.*).cell"),
           (r"fastq.dir/\1.fastq.1.gz",
            r"fastq.dir/\1.fastq.2.gz"))
def cram2fastq(infile, outfiles):
    '''Convert Sanger CRAM files to Fastq format
       Takes care of merging, quality trimming
       and pair reconciliation.
       Intermediate files are not kept by default.'''

    # TODO: make quality trimming optional.

    ###################################
    # set variables and open a log file
    ###################################

    cell_name = os.path.basename(infile)[:-len(".cell")]
    out_dir = os.path.dirname(outfiles[0])
    temp_dir = "fastq.temp.dir"

    log_file = os.path.join(temp_dir,
                            cell_name + ".fastq.extraction.log")

    log = open(log_file, "w")
    log.write("Fastq extraction log file for %(infile)s\n\n")

    def _merge_dicts(a, b):
        x = a.copy()
        x.update(b)
        return(x)

    temp_files = []

    # ##############################################
    # Extract per-end Fastq(s) from the cram file(s)
    # ##############################################

    raw_fastq_names = []
    with open(infile, "rb") as cram_files:

        for line in cram_files:

            cram = line.strip()
            cram_basename = os.path.basename(cram)[:-len(".cram")]
            raw_fastq_name = os.path.join(temp_dir, cram_basename)
            raw_fastq_names.append(raw_fastq_name)

            job_memory = PARAMS["preprocess_memory"]
            statement = '''cramtools fastq
                                       --enumerate
                                       --reverse
                                       -F %(raw_fastq_name)s
                                       -I %(cram)s
                                       --gzip
                         '''
            log.write("Extracting fastqs from %(cram)s:" % locals() + "\n")
            log.write(statement % locals() + "\n")

            P.run(statement)

            log.write("done.\n\n")

    # ####################################
    # Perform quality trimming
    # Merging is also taken care of here.
    # ####################################

    quality = PARAMS["preprocess_quality_threshold"]
    minlen = PARAMS["preprocess_min_length"]
    trim = PARAMS["preprocess_trim"]

    trimmed_fastq_prefix = os.path.join(temp_dir, cell_name)

    trimmed_fastq_files = []

    # fastq(s) for each end are quality trimmed separately
    for end in ["_1", "_2"]:

        raw_fastqs = [x + end + ".fastq.gz" for x in raw_fastq_names]
        temp_files += raw_fastqs

        fastq_list = " ".join(raw_fastqs)

        trimmed_fastq_name = trimmed_fastq_prefix + end + ".trimmed.fastq.gz"
        trimmed_fastq_files.append(trimmed_fastq_name)

        log.write(">> Quality trimming %(fastq_list)s: " % locals() + "\n")
        if trim:
            statement = '''zcat %(fastq_list)s
                       | fastq_quality_trimmer
                           -Q33
                           -t %(quality)s
                           -l %(minlen)s
                       | gzip -c
                       > %(trimmed_fastq_name)s
                       '''
        else:
            statement = '''zcat %(fastq_list)s
                       | gzip -c
                       > %(trimmed_fastq_name)s
                       '''
        log.write(statement % _merge_dicts(PARAMS, locals()) + "\n")
        P.run(statement)
        log.write("done. \n\n")

    # ##################
    # Reconcile the ends
    # ##################

    if PARAMS["preprocess_reconcile"] != "False":

        temp_files += trimmed_fastq_files

        end1, end2 = trimmed_fastq_files

        reconciled_fastq_prefix = outfiles[0][:-len(".1.gz")]

        log.write(">> Reconciling pairs, %(end1)s & %(end2)s: "
                  % locals() + "\n")
        statement = '''python %(scriptsdir)s/fastqs2fastqs.py
                       %(end1)s %(end2)s
                       --method reconcile
                       --chop
                       --unpaired
                       -o "%(reconciled_fastq_prefix)s.%%s.gz";
                    '''
        log.write(statement % _merge_dicts(PARAMS, locals()) + "\n")
        P.run(statement)
        log.write("done\n\n")

    else:
        trimmed_fastq_prefix = outfiles[0][:-len(".1.gz")]
        for end in trimmed_fastq_files:
            if "1.trimmed" in end:
                endn = "1"
            else:
                endn = "2"

            trimmed_end_name = ".".join([trimmed_fastq_prefix, endn, "gz"])
            os.symlink(os.path.abspath(end), trimmed_end_name)

    ##############################
    # Clean up the temporary files
    ##############################

    if PARAMS["keep_temporary"] == 0:

        temp_file_list = " ".join(temp_files)

        # record files sizes and md5 checksums of the temporary files
        log.write(">> Recording sizes and checksums of temporary files:\n")
        statement = '''ls -l %(temp_file_list)s
                       > %(temp_dir)s/%(cell_name)s.ls;
                       checkpoint;
                       md5sum %(temp_file_list)s
                       > %(temp_dir)s/%(cell_name)s.md5;
                    '''
        log.write(statement % locals() + "\n")
        P.run(statement)
        log.write("done\n\n")

        # unlink (delete) the temporary files
        log.write(">> unlinking temporary files: " + temp_file_list + "\n")

        for temp_file in temp_files:
            os.unlink(temp_file)

        log.write("tempororay files unlinked\n")

    log.close()


# ---------------------< generic pipeline tasks >---------------------------- #

@follows(cram2fastq, loadCramQuality)
def full():
    pass


# ########################################################################### #

if __name__ == "__main__":
    sys.exit(P.main(sys.argv))

############################################################################## # #   Kennedy Institute of Rheumatology # #   $Id$ # #   Copyright (C) 2015 <NAME> # #   This program is free software; you can redistribute it and/or #   modify it under the terms of the GNU General Public License #   as published by the Free Software Foundation; either version 2 #   of the License, or (at your option) any later version. # #   This program is distributed in the hope that it will be useful, #   but WITHOUT ANY WARRANTY; without even the implied warranty of #   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the #   GNU General Public License for more details. # #   You should have received a copy of the GNU General Public License #   along with this program; if not, write to the Free Software #   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. ############################################################################### =========================== Pipeline cram2fastq ===========================  :Author: <NAME> :Release: $Id$ :Date: |today| :Tags: Python  Overview ========  This pipeline coverts Sanger CRAM files to fastq.gz, optionally quality trimming and reconciling the fastq files  Usage =====  See :ref:`PipelineSettingUp` and :ref:`PipelineRunning` on general information how to use CGAT pipelines.  Configuration -------------  The pipeline requires a configured :file:`pipeline.ini` file. CGATReport report requires a :file:`conf.py` and optionally a :file:`cgatreport.ini` file (see :ref:`PipelineReporting`).  Default configuration files can be generated by executing:     python <srcdir>/pipeline_cram2fastq.py config  Input files -----------  Requirements ------------   On top of the default CGAT setup, the pipeline requires the following software to be in the path:  .. Add any additional external requirements such as 3rd party software    or R modules below:  Requirements:  * samtools >= 1.1  Pipeline output ===============   Glossary ========   Code ==== # -------------------------- < parse parameters > --------------------------- # # load options from the config file # ----------------------- < pipeline configuration > ------------------------ # # ------------------------< specific pipeline tasks >------------------------ # Validate CRAM files by exit status of        cramtools qstat. Save the quality scores of cram files. temp_quality=`mktemp -p %(cluster_tmpdir)s`;                    cramtools qstat -I %(infile)s > $temp_quality;                    echo $? > %(outfile)s;                    cat $temp_quality                    | awk '{OFS="\\t"} {print $1,$2}'                    > %(outfile_quality)s;                    rm $temp_quality; Check that all crams pass validation or        raise an Error. Load the quality scores for the different cells         into the database (summarized table). Make a table of cells and corresponding cram files # build a dictionary of cell to cram file mappings # write out a per-cell list of cram files Make a per-cell file containing the cram file(s)        corresponding to the cell Convert Sanger CRAM files to Fastq format        Takes care of merging, quality trimming        and pair reconciliation.        Intermediate files are not kept by default. # TODO: make quality trimming optional. ################################### # set variables and open a log file ################################### # ############################################## # Extract per-end Fastq(s) from the cram file(s) # ############################################## cramtools fastq                                        --enumerate                                        --reverse                                        -F %(raw_fastq_name)s                                        -I %(cram)s                                        --gzip # #################################### # Perform quality trimming # Merging is also taken care of here. # #################################### # fastq(s) for each end are quality trimmed separately zcat %(fastq_list)s                        | fastq_quality_trimmer                            -Q33                            -t %(quality)s                            -l %(minlen)s                        | gzip -c                        > %(trimmed_fastq_name)s zcat %(fastq_list)s                        | gzip -c                        > %(trimmed_fastq_name)s # ################## # Reconcile the ends # ################## python %(scriptsdir)s/fastqs2fastqs.py                        %(end1)s %(end2)s                        --method reconcile                        --chop                        --unpaired                        -o "%(reconciled_fastq_prefix)s.%%s.gz"; ############################## # Clean up the temporary files ############################## # record files sizes and md5 checksums of the temporary files ls -l %(temp_file_list)s                        > %(temp_dir)s/%(cell_name)s.ls;                        checkpoint;                        md5sum %(temp_file_list)s                        > %(temp_dir)s/%(cell_name)s.md5; # unlink (delete) the temporary files # ---------------------< generic pipeline tasks >---------------------------- # # ########################################################################### #