<commit_before>"""Provide trimming of input reads from Fastq or BAM files.
"""
import os
import sys

from bcbio.utils import (file_exists, append_stem, replace_directory, is_gzipped)
from bcbio.log import logger
from bcbio.provenance import do
from Bio.Seq import Seq
from bcbio.distributed.transaction import file_transaction

SUPPORTED_ADAPTERS = {
    "illumina": ["AACACTCTTTCCCT", "AGATCGGAAGAGCG"],
    "truseq": ["AGATCGGAAGAG"],
    "polya": ["AAAAAAAAAAAAA"],
    "nextera": ["AATGATACGGCGA", "CAAGCAGAAGACG"]}

def trim_adapters(fastq_files, out_dir, lane_config):
    """
    for small insert sizes, the read length can be longer than the insert
    resulting in the reverse complement of the 3' adapter being sequenced.
    this takes adapter sequences and trims the only the reverse complement
    of the adapter

    MYSEQUENCEAAAARETPADA -> MYSEQUENCEAAAA (no polyA trim)

    """
    quality_format = _get_quality_format(lane_config)
    to_trim = _get_sequences_to_trim(lane_config, SUPPORTED_ADAPTERS)
    out_files = replace_directory(append_stem(fastq_files, ".trimmed"), out_dir)
    out_files = _cutadapt_trim(fastq_files, quality_format, to_trim, out_files)
    return out_files

def _cutadapt_trim(fastq_files, quality_format, adapters, out_files):
    """Trimming with cutadapt, using version installed with bcbio-nextgen.

    Uses the system executable to find the version next to our Anaconda Python.
    TODO: Could we use cutadapt as a library to avoid this?
    """
    if all([file_exists(x) for x in out_files]):
        return out_files
    cmd = _cutadapt_trim_cmd(fastq_files, quality_format, adapters, out_files)
    if len(fastq_files) == 1:
        of1 = out_files[0]
        message = "Trimming %s in single end mode with cutadapt." % (fastq_files[0])
        with file_transaction(of1) as of1_tx:
            do.run(cmd.format(**locals()), message)
    else:
        with file_transaction(out_files) as tx_out_files:
            of1_tx, of2_tx = out_files
            singles_file = of1_tx + ".single"
            message = "Trimming %s and %s in single end mode with cutadapt." % (fastq_files[0],
                                                                                fastq_files[1])
            do.run(cmd.format(**locals()), message)
    return out_files

def _cutadapt_trim_cmd(fastq_files, quality_format, adapters, out_files):
    """Trimming with cutadapt, using version installed with bcbio-nextgen.

    Uses the system executable to find the version next to our Anaconda Python.
    TODO: Could we use cutadapt as a library to avoid this?
    """
    if all([file_exists(x) for x in out_files]):
        return out_files
    if quality_format == "illumina":
        quality_base = "64"
    else:
        quality_base = "33"

    # --times=2 tries twice remove adapters which will allow things like:
    # realsequenceAAAAAAadapter to remove both the poly-A and the adapter
    # this behavior might not be what we want; we could also do two or
    # more passes of cutadapt
    cutadapt = os.path.join(os.path.dirname(sys.executable), "cutadapt")
    adapter_cmd = " ".join(map(lambda x: "--adapter=" + x, adapters))
    base_cmd = ("cutadapt --times=2 --quality-base={quality_base} "
                "--quality-cutoff=5 --format=fastq "
                "--minimum-length=0 {adapter_cmd} ").format(**locals())
    if len(fastq_files) == 1:
        return _cutadapt_se_cmd(fastq_files, out_files, base_cmd)
    else:
        return _cutadapt_pe_cmd(fastq_files, out_files, quality_format, base_cmd)

def _cutadapt_se_cmd(fastq_files, out_files, base_cmd):
    """
    this has to use the -o option, not redirect to stdout in order for gzipping to be
    honored
    """
    fq1 = fastq_files[0]
    of1 = out_files[0]
    cmd = base_cmd + " -o {of1} " + str(fq1)
    return cmd

def _cutadapt_pe_cmd(fastq_files, out_files, quality_format, base_cmd):
    """
    cutadapt can't handle paired end data but sickle can; trim the adapters
    off with cutadapt and then use sickle to trim the reads by length
    """
    qformat = "sanger" if quality_format == "standard" else quality_format
    fq1, fq2 = fastq_files
    first_cmd = "<({base_cmd} {fq1})".format(**locals())
    second_cmd = "<({base_cmd} {fq2})".format(**locals())
    sickle_cmd = ("sickle pe -f {first_cmd} -r {second_cmd} -l 25 -q 0 "
                  "-t {qformat} ")
    sickle_cmd = sickle_cmd + " -g " if is_gzipped(fq1) else sickle_cmd
    sickle_cmd = sickle_cmd.format(**locals())
    sickle_cmd += " -o {of1_tx} -p {of2_tx} -s {singles_file} "
    return sickle_cmd

def _get_sequences_to_trim(lane_config, builtin):
    builtin_adapters = _get_builtin_adapters(lane_config, builtin)
    polya = builtin_adapters.get("polya", [None])[0]
    # allow for trimming of custom sequences for advanced users
    custom_trim = lane_config["algorithm"].get("custom_trim", [])
    builtin_adapters = {k: v for k, v in builtin_adapters.items() if<commit_msg>Remove too-small reads for single-end RNA-seq experiments.<commit_after>"""Provide trimming of input reads from Fastq or BAM files.
"""
import os
import sys

from bcbio.utils import (file_exists, append_stem, replace_directory, is_gzipped)
from bcbio.log import logger
from bcbio.provenance import do
from Bio.Seq import Seq
from bcbio.distributed.transaction import file_transaction

MINIMUM_LENGTH = 25

SUPPORTED_ADAPTERS = {
    "illumina": ["AACACTCTTTCCCT", "AGATCGGAAGAGCG"],
    "truseq": ["AGATCGGAAGAG"],
    "polya": ["AAAAAAAAAAAAA"],
    "nextera": ["AATGATACGGCGA", "CAAGCAGAAGACG"]}

def trim_adapters(fastq_files, out_dir, lane_config):
    """
    for small insert sizes, the read length can be longer than the insert
    resulting in the reverse complement of the 3' adapter being sequenced.
    this takes adapter sequences and trims the only the reverse complement
    of the adapter

    MYSEQUENCEAAAARETPADA -> MYSEQUENCEAAAA (no polyA trim)

    """
    quality_format = _get_quality_format(lane_config)
    to_trim = _get_sequences_to_trim(lane_config, SUPPORTED_ADAPTERS)
    out_files = replace_directory(append_stem(fastq_files, ".trimmed"), out_dir)
    out_files = _cutadapt_trim(fastq_files, quality_format, to_trim, out_files)
    return out_files

def _cutadapt_trim(fastq_files, quality_format, adapters, out_files):
    """Trimming with cutadapt, using version installed with bcbio-nextgen.

    Uses the system executable to find the version next to our Anaconda Python.
    TODO: Could we use cutadapt as a library to avoid this?
    """
    if all([file_exists(x) for x in out_files]):
        return out_files
    cmd = _cutadapt_trim_cmd(fastq_files, quality_format, adapters, out_files)
    if len(fastq_files) == 1:
        of1 = out_files[0]
        message = "Trimming %s in single end mode with cutadapt." % (fastq_files[0])
        with file_transaction(of1) as of1_tx:
            do.run(cmd.format(**locals()), message)
    else:
        with file_transaction(out_files) as tx_out_files:
            of1_tx, of2_tx = out_files
            singles_file = of1_tx + ".single"
            message = "Trimming %s and %s in single end mode with cutadapt." % (fastq_files[0],
                                                                                fastq_files[1])
            do.run(cmd.format(**locals()), message)
    return out_files

def _cutadapt_trim_cmd(fastq_files, quality_format, adapters, out_files):
    """Trimming with cutadapt, using version installed with bcbio-nextgen.

    Uses the system executable to find the version next to our Anaconda Python.
    TODO: Could we use cutadapt as a library to avoid this?
    """
    if all([file_exists(x) for x in out_files]):
        return out_files
    if quality_format == "illumina":
        quality_base = "64"
    else:
        quality_base = "33"

    # --times=2 tries twice remove adapters which will allow things like:
    # realsequenceAAAAAAadapter to remove both the poly-A and the adapter
    # this behavior might not be what we want; we could also do two or
    # more passes of cutadapt
    cutadapt = os.path.join(os.path.dirname(sys.executable), "cutadapt")
    adapter_cmd = " ".join(map(lambda x: "--adapter=" + x, adapters))
    base_cmd = ("cutadapt --times=2 --quality-base={quality_base} "
                "--quality-cutoff=5 --format=fastq "
                "{adapter_cmd} ").format(**locals())
    if len(fastq_files) == 1:
        return _cutadapt_se_cmd(fastq_files, out_files, base_cmd)
    else:
        return _cutadapt_pe_cmd(fastq_files, out_files, quality_format, base_cmd)

def _cutadapt_se_cmd(fastq_files, out_files, base_cmd):
    """
    this has to use the -o option, not redirect to stdout in order for gzipping to be
    honored
    """
    min_length = MINIMUM_LENGTH
    cmd = base_cmd + " --minimum-length={min_length} ".format(**locals())
    fq1 = fastq_files[0]
    of1 = out_files[0]
    cmd += " -o {of1} " + str(fq1)
    return cmd

def _cutadapt_pe_cmd(fastq_files, out_files, quality_format, base_cmd):
    """
    cutadapt can't handle paired end data but sickle can; trim the adapters
    off with cutadapt and then use sickle to trim the reads by length
    """
    min_length = MINIMUM_LENGTH
    qformat = "sanger" if quality_format == "standard" else quality_format
    fq1, fq2 = fastq_files
    base_cmd += " --minimum-length=0 "
    first_cmd = "<({base_cmd} {fq1})".format(**locals())
    second_cmd = "<({base_cmd} {fq2})".format(**locals())
    sickle_cmd = ("sickle pe -f {first_cmd} -r {second_cmd} -l {min_length} -q 0 "
                  "-t {qformat} ")
    sickle_cmd = sickle_cmd + " -g " if is_gzipped(fq1) else sickle_cmd
    sickle_cmd = sickle_cmd.format(**locals())
    sickle_cmd += " -o {of1_tx} -p {of2_tx} -s {singles_file} "
    return sickle_cmd

def _get_sequences_to_trim(lane_config, builtin):
    builtin_adapters = _get_builtin_adapters(lane_config, builtin)
    polya = builtin_adapters.get("polya", [None])[0]
    # allow for trimming of custom sequences for advanced users
    custom_trim = lane_config["algorithm"].get("custom_trim", [])
    builtin_adapters = {k: v for k, v in builtin_adapters.items() if
{'fraction_non_alphanumeric': 0.1014720314033366, 'fraction_numerical': 0.00745829244357213, 'mean_word_length': 3.3182203389830507, 'pattern_counts': {'":': 10, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5957551', 'n_tokens_mistral': 3451, 'n_tokens_neox': 3310, 'n_words': 955}