MikkelSchubert / adapterremoval / #45

Committed 20 Sep 2024 06:49PM UTC coverage: 26.244% (-49.2%) from 75.443%

Build # #45

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

attempt to fix coveralls run

Run Details

2458 of 9366 relevant lines covered (26.24%)

4362.23 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/trimming.cpp

/*************************************************************************\
 * AdapterRemoval - cleaning next-generation sequencing reads            *
 *                                                                       *
 * Copyright (C) 2011 by Stinus Lindgreen - stinus@binf.ku.dk            *
 * Copyright (C) 2014 by Mikkel Schubert - mikkelsch@gmail.com           *
 *                                                                       *
 * 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 3 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, see <http://www.gnu.org/licenses/>. *
\*************************************************************************/
#include "trimming.hpp"
#include "alignment.hpp"   // for alignment_info, sequence_merger, sequence_...
#include "commontypes.hpp" // for read_type, trimming_strategy, merge_strategy
#include "counts.hpp"      // for counts, indexed_count
#include "debug.hpp"       // for AR_FAIL, AR_REQUIRE
#include "fastq_io.hpp"    // for chunk_ptr, fastq_...
#include "output.hpp"      // for sample_output_files, processed_reads
#include "sequence_sets.hpp" // for adapter_set
#include "serializer.hpp"    // for fastq_flags
#include "simd.hpp"          // for size_t
#include "statistics.hpp" // for trimming_statistics, reads_and_bases, fast...
#include "userconfig.hpp" // for userconfig
#include <cstddef>        // for size_t
#include <memory>         // for unique_ptr, __shared_ptr_access, make_unique
#include <string>         // for string
#include <utility>        // for pair, move

namespace adapterremoval {

////////////////////////////////////////////////////////////////////////////////
// Helper functions

namespace {

/** Tracks overlapping reads, to account for trimming affecting the overlap */
class merged_reads
{
public:
  merged_reads(const fastq& read1, const fastq& read2, size_t insert_size)
    : m_len_1(read1.length())
    , m_len_2(read2.length())
    , m_overlap(m_len_1 + m_len_2 - insert_size)
  {
  }

  /** Increments bases trimmed and returns true if overlap was trimmed */
  bool increment(const size_t trim5p, const size_t trim3p)
  {
    if (m_trimmed_5p + m_trimmed_3p < m_len_1 + m_len_2 - m_overlap) {
      m_trimmed_5p += trim5p;
      m_trimmed_3p += trim3p;

      return (trim5p && trim3p) ||
             ((trim5p || trim3p) && ((m_trimmed_5p > m_len_1 - m_overlap) ||
                                     (m_trimmed_3p > m_len_2 - m_overlap)));
    } else {
      return false;
    }
  }

private:
  const size_t m_len_1;
  const size_t m_len_2;
  const size_t m_overlap;

  size_t m_trimmed_5p = 0;
  size_t m_trimmed_3p = 0;
};

/** Trims poly-X tails from sequence prior to adapter trimming **/
void
pre_trim_poly_x_tail(const userconfig& config,
                     trimming_statistics& stats,
                     fastq& read)
{
  if (!config.pre_trim_poly_x.empty()) {
    const auto result = read.poly_x_trimming(config.pre_trim_poly_x,
                                             config.trim_poly_x_threshold);

    if (result.second) {
      stats.poly_x_pre_trimmed_reads.inc(result.first);
      stats.poly_x_pre_trimmed_bases.inc(result.first, result.second);
    }
  }
}

/** Trims poly-X tails from sequence after adapter trimming **/
void
post_trim_poly_x_tail(const userconfig& config,
                      trimming_statistics& stats,
                      fastq& read)
{
  if (!config.post_trim_poly_x.empty()) {
    const auto result = read.poly_x_trimming(config.post_trim_poly_x,
                                             config.trim_poly_x_threshold);

    if (result.second) {
      stats.poly_x_post_trimmed_reads.inc(result.first);
      stats.poly_x_post_trimmed_bases.inc(result.first, result.second);
    }
  }
}

/** Trims fixed bases from read termini and returns the number trimmed. **/
void
trim_read_termini(reads_and_bases& stats,
                  fastq& read,
                  size_t trim_5p,
                  size_t trim_3p)
{
  const auto length = read.length();
  if ((trim_5p || trim_3p) && length) {
    if (trim_5p + trim_3p < length) {
      read.truncate(trim_5p, length - trim_5p - trim_3p);
    } else {
      read.truncate(0, 0);
    }

    stats.inc(length - read.length());
  }
}

/** Trims fixed number of 5'/3' bases prior to adapter trimming */
void
pre_trim_read_termini(const userconfig& config,
                      trimming_statistics& stats,
                      fastq& read,
                      read_type type)
{
  size_t trim_5p = 0;
  size_t trim_3p = 0;

  if (type == read_type::mate_1) {
#ifdef PRE_TRIM_5P
    trim_5p = config.pre_trim_fixed_5p.first;
#endif
    trim_3p = config.pre_trim_fixed_3p.first;
  } else if (type == read_type::mate_2) {
#ifdef PRE_TRIM_5P
    trim_5p = config.pre_trim_fixed_5p.second;
#endif
    trim_3p = config.pre_trim_fixed_3p.second;
  } else {
    AR_FAIL("invalid read type in pre_trim_read_termini");
  }

  trim_read_termini(stats.terminal_pre_trimmed, read, trim_5p, trim_3p);
}

/** Trims fixed number of 5'/3' bases after adapter trimming */
void
post_trim_read_termini(const userconfig& config,
                       trimming_statistics& stats,
                       fastq& read,
                       read_type type,
                       merged_reads* mstats = nullptr)
{
  size_t trim_5p = 0;
  size_t trim_3p = 0;

  switch (type) {
    case read_type::mate_1:
      trim_5p = config.post_trim_fixed_5p.first;
      trim_3p = config.post_trim_fixed_3p.first;
      break;

    case read_type::mate_2:
      trim_5p = config.post_trim_fixed_5p.second;
      trim_3p = config.post_trim_fixed_3p.second;
      break;

    case read_type::merged:
      AR_REQUIRE(config.paired_ended_mode);
      trim_5p = config.post_trim_fixed_5p.first;
      trim_3p = config.post_trim_fixed_5p.second;
      break;

    case read_type::singleton:
    case read_type::discarded:
      AR_FAIL("unsupported read type in post_trim_read_termini");

    case read_type::max:
    default:
      AR_FAIL("invalid read type in post_trim_read_termini");
  }

  trim_read_termini(stats.terminal_post_trimmed, read, trim_5p, trim_3p);
  if (mstats && mstats->increment(trim_5p, trim_3p)) {
    stats.terminal_post_trimmed.inc_reads();
  }
}

/** Quality trims a read after adapter trimming */
void
post_trim_read_by_quality(const userconfig& config,
                          trimming_statistics& stats,
                          fastq& read,
                          merged_reads* mstats = nullptr)
{
  fastq::ntrimmed trimmed;
  switch (config.trim) {
    case trimming_strategy::none:
      break;

    case trimming_strategy::mott:
      trimmed = read.mott_trimming(config.trim_mott_rate, config.preserve5p);
      break;

    case trimming_strategy::window:
      trimmed = read.trim_windowed_bases(config.trim_ambiguous_bases,
                                         config.trim_quality_score,
                                         config.trim_window_length,
                                         config.preserve5p);
      break;

    case trimming_strategy::per_base:
      trimmed = read.trim_trailing_bases(
        config.trim_ambiguous_bases,
        config.trim_low_quality_bases ? config.trim_quality_score : -1,
        config.preserve5p);
      break;

    default:
      AR_FAIL("not implemented");
  }

  if (trimmed.first || trimmed.second) {
    stats.low_quality_trimmed.inc(trimmed.first + trimmed.second);

    if (mstats && mstats->increment(trimmed.first, trimmed.second)) {
      stats.low_quality_trimmed.inc_reads();
    }
  }
}

bool
is_acceptable_read(const userconfig& config,
                   trimming_statistics& stats,
                   const fastq& seq,
                   const size_t n_reads = 1)
{
  const auto length = seq.length();

  if (length < config.min_genomic_length) {
    stats.filtered_min_length.inc_reads(n_reads);
    stats.filtered_min_length.inc_bases(length);
    return false;
  } else if (length > config.max_genomic_length) {
    stats.filtered_max_length.inc_reads(n_reads);
    stats.filtered_max_length.inc_bases(length);
    return false;
  }

  const auto max_n = config.max_ambiguous_bases;
  if (max_n < length && seq.count_ns() > max_n) {
    stats.filtered_ambiguous.inc_reads(n_reads);
    stats.filtered_ambiguous.inc_bases(length);
    return false;
  }

  if (length > 0 && config.min_mean_quality > 0.0 &&
      seq.mean_quality() < config.min_mean_quality) {
    stats.filtered_mean_quality.inc_reads(n_reads);
    stats.filtered_mean_quality.inc_bases(length);
    return false;
  }

  if (config.min_complexity > 0.0 && seq.complexity() < config.min_complexity) {
    stats.filtered_low_complexity.inc_reads(n_reads);
    stats.filtered_low_complexity.inc_bases(length);
    return false;
  }

  return true;
}

} // namespace

////////////////////////////////////////////////////////////////////////////////
// Implementations for `reads_processor`

reads_processor::reads_processor(const userconfig& config,
                                 const sample_output_files& output,
                                 const size_t nth,
                                 trim_stats_ptr sink)
  : analytical_step(processing_order::unordered, "reads_processor")
  , m_config(config)
  , m_output(output)
  , m_sample(nth)
  , m_stats_sink(std::move(sink))
{
  AR_REQUIRE(m_stats_sink);

  m_stats.emplace_back_n(m_config.max_threads, m_config.report_sample_rate);
}

void
reads_processor::finalize()
{
  m_stats.merge_into(*m_stats_sink);
}

////////////////////////////////////////////////////////////////////////////////
// Implementations for `se_reads_processor`

se_reads_processor::se_reads_processor(const userconfig& config,
                                       const sample_output_files& output,
                                       const size_t nth,
                                       trim_stats_ptr sink)
  : reads_processor(config, output, nth, sink)
{
}

chunk_vec
se_reads_processor::process(chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  processed_reads chunks{ m_output };
  chunks.set_sample(m_config.samples.at(m_sample));
  chunks.set_mate_separator(chunk->mate_separator);

  if (chunk->first) {
    chunks.write_headers(m_config.args);
  }

  auto stats = m_stats.acquire();
  stats->adapter_trimmed_reads.resize_up_to(m_config.samples.adapters().size());
  stats->adapter_trimmed_bases.resize_up_to(m_config.samples.adapters().size());

  // A sequence aligner per barcode (pair)
  std::vector<sequence_aligner> aligners;
  for (const auto& it : m_config.samples.at(m_sample)) {
    aligners.emplace_back(it.adapters, m_config.simd);
  }

  AR_REQUIRE(!aligners.empty());
  AR_REQUIRE(chunk->barcodes.empty() ||
             chunk->barcodes.size() == chunk->reads_1.size());

  auto barcode_it = chunk->barcodes.begin();
  const auto barcode_end = chunk->barcodes.end();

  for (auto& read : chunk->reads_1) {
    const size_t in_length = read.length();

    // Trim fixed number of bases from 5' and/or 3' termini
    pre_trim_read_termini(m_config, *stats, read, read_type::mate_1);
    // Trim poly-X tails for zero or more X
    pre_trim_poly_x_tail(m_config, *stats, read);

    const auto barcode = (barcode_it == barcode_end) ? 0 : *barcode_it++;
    const auto alignment =
      aligners.at(barcode).align_single_end(read, m_config.shift);

    if (m_config.is_good_alignment(alignment)) {
      const auto length = read.length();
      alignment.truncate_single_end(read);

      stats->adapter_trimmed_reads.inc(alignment.adapter_id);
      stats->adapter_trimmed_bases.inc(alignment.adapter_id,
                                       length - read.length());
    }

    // Add (optional) user specified prefixes to read names
    read.add_prefix_to_name(m_config.prefix_read_1);

    // Trim fixed number of bases from 5' and/or 3' termini
    post_trim_read_termini(m_config, *stats, read, read_type::mate_1);
    // Trim poly-X tails for zero or more X
    post_trim_poly_x_tail(m_config, *stats, read);
    // Sliding window trimming or single-base trimming of low quality bases
    post_trim_read_by_quality(m_config, *stats, read);

    stats->total_trimmed.inc_reads(in_length != read.length());
    stats->total_trimmed.inc_bases(in_length - read.length());

    if (is_acceptable_read(m_config, *stats, read)) {
      stats->read_1->process(read);
      chunks.add(read, read_type::mate_1, fastq_flags::se, barcode);
    } else {
      stats->discarded->process(read);
      chunks.add(read, read_type::discarded, fastq_flags::se_fail, barcode);
    }
  }

  m_stats.release(stats);

  return chunks.finalize(chunk->eof);
}

////////////////////////////////////////////////////////////////////////////////
// Implementations for `pe_reads_processor`

void
add_pe_statistics(const trimming_statistics& stats,
                  const fastq& read,
                  read_type type)
{
  switch (type) {
    case read_type::mate_1:
      stats.read_1->process(read);
      break;
    case read_type::mate_2:
      stats.read_2->process(read);
      break;

    case read_type::merged:
      stats.merged->process(read);
      break;

    case read_type::singleton:
      stats.singleton->process(read);
      break;

    case read_type::discarded:
      stats.discarded->process(read);
      break;

    case read_type::max:
      AR_FAIL("unhandled read type");

    default:
      AR_FAIL("invalid read type");
  }
}

pe_reads_processor::pe_reads_processor(const userconfig& config,
                                       const sample_output_files& output,
                                       const size_t nth,
                                       trim_stats_ptr sink)
  : reads_processor(config, output, nth, sink)
{
}

chunk_vec
pe_reads_processor::process(chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  processed_reads chunks{ m_output };
  chunks.set_sample(m_config.samples.at(m_sample));
  chunks.set_mate_separator(chunk->mate_separator);

  if (chunk->first) {
    chunks.write_headers(m_config.args);
  }

  sequence_merger merger;
  merger.set_merge_strategy(m_config.merge);
  merger.set_max_recalculated_score(m_config.merge_quality_max);

  // A sequence aligner per barcode (pair)
  std::vector<sequence_aligner> aligners;
  for (const auto& it : m_config.samples.at(m_sample)) {
    aligners.emplace_back(it.adapters, m_config.simd);
  }

  auto stats = m_stats.acquire();
  stats->adapter_trimmed_reads.resize_up_to(m_config.samples.adapters().size());
  stats->adapter_trimmed_bases.resize_up_to(m_config.samples.adapters().size());

  AR_REQUIRE(!aligners.empty());
  AR_REQUIRE(chunk->reads_1.size() == chunk->reads_2.size());
  AR_REQUIRE(chunk->barcodes.empty() ||
             chunk->barcodes.size() == chunk->reads_1.size());

  auto barcode_it = chunk->barcodes.begin();
  const auto barcode_end = chunk->barcodes.end();

  auto it_1 = chunk->reads_1.begin();
  auto it_2 = chunk->reads_2.begin();
  while (it_1 != chunk->reads_1.end()) {
    fastq& read_1 = *it_1++;
    fastq& read_2 = *it_2++;

    const size_t in_length_1 = read_1.length();
    const size_t in_length_2 = read_2.length();

    // Trim fixed number of bases from 5' and/or 3' termini
    pre_trim_read_termini(m_config, *stats, read_1, read_type::mate_1);
    pre_trim_read_termini(m_config, *stats, read_2, read_type::mate_2);

    // Trim poly-X tails for zero or more X
    pre_trim_poly_x_tail(m_config, *stats, read_1);
    pre_trim_poly_x_tail(m_config, *stats, read_2);

    // Reverse complement to match the orientation of read_1
    read_2.reverse_complement();

    const auto barcode = (barcode_it == barcode_end) ? 0 : *barcode_it++;
    const auto alignment =
      aligners.at(barcode).align_paired_end(read_1, read_2, m_config.shift);

    if (m_config.is_good_alignment(alignment)) {
      const size_t insert_size = alignment.insert_size(read_1, read_2);
      const bool can_merge_alignment = m_config.can_merge_alignment(alignment);
      if (can_merge_alignment) {
        // Insert size calculated from untrimmed reads
        stats->insert_sizes.resize_up_to(insert_size + 1);
        stats->insert_sizes.inc(insert_size);
      }

      const size_t pre_trimmed_bp = read_1.length() + read_2.length();
      const size_t n_adapters = alignment.truncate_paired_end(read_1, read_2);
      const size_t post_trimmed_bp = read_1.length() + read_2.length();

      stats->adapter_trimmed_reads.inc(alignment.adapter_id, n_adapters);
      stats->adapter_trimmed_bases.inc(alignment.adapter_id,
                                       pre_trimmed_bp - post_trimmed_bp);

      if (m_config.merge != merge_strategy::none && can_merge_alignment) {
        // Track if one or both source reads are trimmed post merging
        merged_reads mstats(read_1, read_2, insert_size);

        // Merge read_2 into read_1
        merger.merge(alignment, read_1, read_2);

        // Add (optional) user specified prefix to read names
        read_1.add_prefix_to_name(m_config.prefix_merged);

        // Trim fixed number of bases from 5' and/or 3' termini
        post_trim_read_termini(
          m_config, *stats, read_1, read_type::merged, &mstats);

        if (!m_config.preserve5p) {
          // A merged read essentially consists of two 5p termini, so neither
          // end should be trimmed if --preserve5p is used.
          post_trim_read_by_quality(m_config, *stats, read_1, &mstats);
        }

        // Track total amount of bases trimmed/lost
        stats->total_trimmed.inc_reads(2);
        stats->total_trimmed.inc_bases(in_length_1 + in_length_2 -
                                       read_1.length());

        if (is_acceptable_read(m_config, *stats, read_1, 2)) {
          stats->merged->process(read_1, 2);
          chunks.add(read_1, read_type::merged, fastq_flags::se, barcode);
        } else {
          stats->discarded->process(read_1, 2);
          chunks.add(
            read_1, read_type::discarded, fastq_flags::se_fail, barcode);
        }

        continue;
      }
    }

    // Reads were not aligned or merging is not enabled
    // Undo reverse complementation (post truncation of adapters)
    read_2.reverse_complement();

    // Add (optional) user specified prefixes to read names
    read_1.add_prefix_to_name(m_config.prefix_read_1);
    read_2.add_prefix_to_name(m_config.prefix_read_2);

    // Trim fixed number of bases from 5' and/or 3' termini
    post_trim_read_termini(m_config, *stats, read_1, read_type::mate_1);
    post_trim_read_termini(m_config, *stats, read_2, read_type::mate_2);

    // Trim poly-X tails for zero or more X
    post_trim_poly_x_tail(m_config, *stats, read_1);
    post_trim_poly_x_tail(m_config, *stats, read_2);

    // Sliding window trimming or single-base trimming of low quality bases
    post_trim_read_by_quality(m_config, *stats, read_1);
    post_trim_read_by_quality(m_config, *stats, read_2);

    // Are the reads good enough? Not too many Ns?
    const bool is_ok_1 = is_acceptable_read(m_config, *stats, read_1);
    const bool is_ok_2 = is_acceptable_read(m_config, *stats, read_2);

    read_type type_1;
    read_type type_2;
    if (is_ok_1 && is_ok_2) {
      type_1 = read_type::mate_1;
      type_2 = read_type::mate_2;
    } else if (is_ok_1) {
      type_1 = read_type::singleton;
      type_2 = read_type::discarded;
    } else if (is_ok_2) {
      type_1 = read_type::discarded;
      type_2 = read_type::singleton;
    } else {
      type_1 = read_type::discarded;
      type_2 = read_type::discarded;
    }

    add_pe_statistics(*stats, read_1, type_1);
    add_pe_statistics(*stats, read_2, type_2);

    // Total number of reads/bases trimmed
    stats->total_trimmed.inc_reads((in_length_1 != read_1.length()) +
                                   (in_length_2 != read_2.length()));
    stats->total_trimmed.inc_bases((in_length_1 - read_1.length()) +
                                   (in_length_2 - read_2.length()));

    const auto flags_1 = is_ok_1 ? fastq_flags::pe_1 : fastq_flags::pe_1_fail;
    const auto flags_2 = is_ok_2 ? fastq_flags::pe_2 : fastq_flags::pe_2_fail;

    // Queue reads last, since this result in modifications to lengths
    chunks.add(read_1, type_1, flags_1, barcode);
    chunks.add(read_2, type_2, flags_2, barcode);
  }

  // Track amount of overlapping bases "lost" due to read merging
  stats->reads_merged.inc_reads(merger.reads_merged());
  stats->reads_merged.inc_bases(merger.bases_merged());

  m_stats.release(stats);

  return chunks.finalize(chunk->eof);
}

} // namespace adapterremoval

1	/*************************************************************************\
2	* AdapterRemoval - cleaning next-generation sequencing reads *
3	* *
4	* Copyright (C) 2011 by Stinus Lindgreen - stinus@binf.ku.dk *
5	* Copyright (C) 2014 by Mikkel Schubert - mikkelsch@gmail.com *
6	* *
7	* This program is free software: you can redistribute it and/or modify *
8	* it under the terms of the GNU General Public License as published by *
9	* the Free Software Foundation, either version 3 of the License, or *
10	* (at your option) any later version. *
11	* *
12	* This program is distributed in the hope that it will be useful, *
13	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15	* GNU General Public License for more details. *
16	* *
17	* You should have received a copy of the GNU General Public License *
18	* along with this program. If not, see <http://www.gnu.org/licenses/>. *
19	\*************************************************************************/
20	#include "trimming.hpp"
21	#include "alignment.hpp" // for alignment_info, sequence_merger, sequence_...
22	#include "commontypes.hpp" // for read_type, trimming_strategy, merge_strategy
23	#include "counts.hpp" // for counts, indexed_count
24	#include "debug.hpp" // for AR_FAIL, AR_REQUIRE
25	#include "fastq_io.hpp" // for chunk_ptr, fastq_...
26	#include "output.hpp" // for sample_output_files, processed_reads
27	#include "sequence_sets.hpp" // for adapter_set
28	#include "serializer.hpp" // for fastq_flags
29	#include "simd.hpp" // for size_t
30	#include "statistics.hpp" // for trimming_statistics, reads_and_bases, fast...
31	#include "userconfig.hpp" // for userconfig
32	#include <cstddef> // for size_t
33	#include <memory> // for unique_ptr, __shared_ptr_access, make_unique
34	#include <string> // for string
35	#include <utility> // for pair, move
36
37	namespace adapterremoval {
38
39	////////////////////////////////////////////////////////////////////////////////
40	// Helper functions
41
42	namespace {
43
44	/** Tracks overlapping reads, to account for trimming affecting the overlap */
45	class merged_reads
46	{
47	public:
48	merged_reads(const fastq& read1, const fastq& read2, size_t insert_size)	×
49	: m_len_1(read1.length())	×
50	, m_len_2(read2.length())	×
51	, m_overlap(m_len_1 + m_len_2 - insert_size)	×
52	{
53	}
54
55	/** Increments bases trimmed and returns true if overlap was trimmed */
56	bool increment(const size_t trim5p, const size_t trim3p)	×
57	{
58	if (m_trimmed_5p + m_trimmed_3p < m_len_1 + m_len_2 - m_overlap) {	×
59	m_trimmed_5p += trim5p;	×
60	m_trimmed_3p += trim3p;	×
61
62	return (trim5p && trim3p) \|\|	×
63	((trim5p \|\| trim3p) && ((m_trimmed_5p > m_len_1 - m_overlap) \|\|	×
64	(m_trimmed_3p > m_len_2 - m_overlap)));	×
65	} else {
66	return false;
67	}
68	}
69
70	private:
71	const size_t m_len_1;
72	const size_t m_len_2;
73	const size_t m_overlap;
74
75	size_t m_trimmed_5p = 0;
76	size_t m_trimmed_3p = 0;
77	};
78
79	/ Trims poly-X tails from sequence prior to adapter trimming /
80	void
81	pre_trim_poly_x_tail(const userconfig& config,	×
82	trimming_statistics& stats,
83	fastq& read)
84	{
85	if (!config.pre_trim_poly_x.empty()) {	×
86	const auto result = read.poly_x_trimming(config.pre_trim_poly_x,	×
87	config.trim_poly_x_threshold);	×
88
89	if (result.second) {	×
90	stats.poly_x_pre_trimmed_reads.inc(result.first);	×
91	stats.poly_x_pre_trimmed_bases.inc(result.first, result.second);	×
92	}
93	}
94	}
95
96	/ Trims poly-X tails from sequence after adapter trimming /
97	void
98	post_trim_poly_x_tail(const userconfig& config,	×
99	trimming_statistics& stats,
100	fastq& read)
101	{
102	if (!config.post_trim_poly_x.empty()) {	×
103	const auto result = read.poly_x_trimming(config.post_trim_poly_x,	×
104	config.trim_poly_x_threshold);	×
105
106	if (result.second) {	×
107	stats.poly_x_post_trimmed_reads.inc(result.first);	×
108	stats.poly_x_post_trimmed_bases.inc(result.first, result.second);	×
109	}
110	}
111	}
112
113	/ Trims fixed bases from read termini and returns the number trimmed. /
114	void
115	trim_read_termini(reads_and_bases& stats,	×
116	fastq& read,
117	size_t trim_5p,
118	size_t trim_3p)
119	{
120	const auto length = read.length();	×
121	if ((trim_5p \|\| trim_3p) && length) {	×
122	if (trim_5p + trim_3p < length) {	×
123	read.truncate(trim_5p, length - trim_5p - trim_3p);	×
124	} else {
125	read.truncate(0, 0);	×
126	}
127
128	stats.inc(length - read.length());	×
129	}
130	}
131
132	/** Trims fixed number of 5'/3' bases prior to adapter trimming */
133	void
134	pre_trim_read_termini(const userconfig& config,	×
135	trimming_statistics& stats,
136	fastq& read,
137	read_type type)
138	{
139	size_t trim_5p = 0;	×
140	size_t trim_3p = 0;	×
141
142	if (type == read_type::mate_1) {	×
143	#ifdef PRE_TRIM_5P
144	trim_5p = config.pre_trim_fixed_5p.first;
145	#endif
146	trim_3p = config.pre_trim_fixed_3p.first;	×
147	} else if (type == read_type::mate_2) {	×
148	#ifdef PRE_TRIM_5P
149	trim_5p = config.pre_trim_fixed_5p.second;
150	#endif
151	trim_3p = config.pre_trim_fixed_3p.second;	×
152	} else {
153	AR_FAIL("invalid read type in pre_trim_read_termini");
154	}
155
156	trim_read_termini(stats.terminal_pre_trimmed, read, trim_5p, trim_3p);	×
157	}
158
159	/** Trims fixed number of 5'/3' bases after adapter trimming */
160	void
161	post_trim_read_termini(const userconfig& config,	×
162	trimming_statistics& stats,
163	fastq& read,
164	read_type type,
165	merged_reads* mstats = nullptr)
166	{
167	size_t trim_5p = 0;	×
168	size_t trim_3p = 0;	×
169
170	switch (type) {	×
171	case read_type::mate_1:	×
172	trim_5p = config.post_trim_fixed_5p.first;	×
173	trim_3p = config.post_trim_fixed_3p.first;	×
174	break;	×
175
176	case read_type::mate_2:	×
177	trim_5p = config.post_trim_fixed_5p.second;	×
178	trim_3p = config.post_trim_fixed_3p.second;	×
179	break;	×
180
181	case read_type::merged:	×
182	AR_REQUIRE(config.paired_ended_mode);	×
183	trim_5p = config.post_trim_fixed_5p.first;	×
184	trim_3p = config.post_trim_fixed_5p.second;	×
185	break;	×
186
187	case read_type::singleton:	×
188	case read_type::discarded:	×
189	AR_FAIL("unsupported read type in post_trim_read_termini");	×
190
191	case read_type::max:	×
192	default:	×
193	AR_FAIL("invalid read type in post_trim_read_termini");	×
194	}
195
196	trim_read_termini(stats.terminal_post_trimmed, read, trim_5p, trim_3p);	×
197	if (mstats && mstats->increment(trim_5p, trim_3p)) {	×
198	stats.terminal_post_trimmed.inc_reads();	×
199	}
200	}
201
202	/** Quality trims a read after adapter trimming */
203	void
204	post_trim_read_by_quality(const userconfig& config,	×
205	trimming_statistics& stats,
206	fastq& read,
207	merged_reads* mstats = nullptr)
208	{
209	fastq::ntrimmed trimmed;	×
210	switch (config.trim) {	×
211	case trimming_strategy::none:
212	break;
213
214	case trimming_strategy::mott:	×
215	trimmed = read.mott_trimming(config.trim_mott_rate, config.preserve5p);	×
216	break;	×
217
218	case trimming_strategy::window:	×
219	trimmed = read.trim_windowed_bases(config.trim_ambiguous_bases,	×
220	config.trim_quality_score,	×
221	config.trim_window_length,	×
222	config.preserve5p);	×
223	break;	×
224
225	case trimming_strategy::per_base:	×
226	trimmed = read.trim_trailing_bases(	×
227	config.trim_ambiguous_bases,	×
228	config.trim_low_quality_bases ? config.trim_quality_score : -1,	×
229	config.preserve5p);	×
230	break;	×
231
232	default:	×
233	AR_FAIL("not implemented");	×
234	}
235
236	if (trimmed.first \|\| trimmed.second) {	×
237	stats.low_quality_trimmed.inc(trimmed.first + trimmed.second);	×
238
239	if (mstats && mstats->increment(trimmed.first, trimmed.second)) {	×
240	stats.low_quality_trimmed.inc_reads();	×
241	}
242	}
243	}
244
245	bool
246	is_acceptable_read(const userconfig& config,	×
247	trimming_statistics& stats,
248	const fastq& seq,
249	const size_t n_reads = 1)
250	{
251	const auto length = seq.length();	×
252
253	if (length < config.min_genomic_length) {	×
254	stats.filtered_min_length.inc_reads(n_reads);	×
255	stats.filtered_min_length.inc_bases(length);	×
256	return false;	×
257	} else if (length > config.max_genomic_length) {	×
258	stats.filtered_max_length.inc_reads(n_reads);	×
259	stats.filtered_max_length.inc_bases(length);	×
260	return false;	×
261	}
262
263	const auto max_n = config.max_ambiguous_bases;	×
264	if (max_n < length && seq.count_ns() > max_n) {	×
265	stats.filtered_ambiguous.inc_reads(n_reads);	×
266	stats.filtered_ambiguous.inc_bases(length);	×
267	return false;	×
268	}
269
270	if (length > 0 && config.min_mean_quality > 0.0 &&	×
271	seq.mean_quality() < config.min_mean_quality) {	×
272	stats.filtered_mean_quality.inc_reads(n_reads);	×
273	stats.filtered_mean_quality.inc_bases(length);	×
274	return false;	×
275	}
276
277	if (config.min_complexity > 0.0 && seq.complexity() < config.min_complexity) {	×
278	stats.filtered_low_complexity.inc_reads(n_reads);	×
279	stats.filtered_low_complexity.inc_bases(length);	×
280	return false;	×
281	}
282
283	return true;
284	}
285
286	} // namespace
287
288	////////////////////////////////////////////////////////////////////////////////
289	// Implementations for `reads_processor`
290
291	reads_processor::reads_processor(const userconfig& config,	×
292	const sample_output_files& output,
293	const size_t nth,
294	trim_stats_ptr sink)	×
295	: analytical_step(processing_order::unordered, "reads_processor")
296	, m_config(config)	×
297	, m_output(output)	×
298	, m_sample(nth)	×
299	, m_stats_sink(std::move(sink))	×
300	{
301	AR_REQUIRE(m_stats_sink);	×
302
303	m_stats.emplace_back_n(m_config.max_threads, m_config.report_sample_rate);	×
304	}
305
306	void
307	reads_processor::finalize()	×
308	{
309	m_stats.merge_into(*m_stats_sink);	×
310	}
311
312	////////////////////////////////////////////////////////////////////////////////
313	// Implementations for `se_reads_processor`
314
315	se_reads_processor::se_reads_processor(const userconfig& config,	×
316	const sample_output_files& output,
317	const size_t nth,
318	trim_stats_ptr sink)	×
319	: reads_processor(config, output, nth, sink)	×
320	{
321	}
322
323	chunk_vec
324	se_reads_processor::process(chunk_ptr chunk)	×
325	{
326	AR_REQUIRE(chunk);	×
327	processed_reads chunks{ m_output };	×
328	chunks.set_sample(m_config.samples.at(m_sample));	×
329	chunks.set_mate_separator(chunk->mate_separator);	×
330
331	if (chunk->first) {	×
332	chunks.write_headers(m_config.args);	×
333	}
334
335	auto stats = m_stats.acquire();	×
336	stats->adapter_trimmed_reads.resize_up_to(m_config.samples.adapters().size());	×
337	stats->adapter_trimmed_bases.resize_up_to(m_config.samples.adapters().size());	×
338
339	// A sequence aligner per barcode (pair)
340	std::vector<sequence_aligner> aligners;	×
341	for (const auto& it : m_config.samples.at(m_sample)) {	×
342	aligners.emplace_back(it.adapters, m_config.simd);	×
343	}
344
345	AR_REQUIRE(!aligners.empty());	×
346	AR_REQUIRE(chunk->barcodes.empty() \|\|	×
347	chunk->barcodes.size() == chunk->reads_1.size());
348
349	auto barcode_it = chunk->barcodes.begin();	×
350	const auto barcode_end = chunk->barcodes.end();	×
351
352	for (auto& read : chunk->reads_1) {	×
353	const size_t in_length = read.length();	×
354
355	// Trim fixed number of bases from 5' and/or 3' termini
356	pre_trim_read_termini(m_config, *stats, read, read_type::mate_1);	×
357	// Trim poly-X tails for zero or more X
358	pre_trim_poly_x_tail(m_config, *stats, read);	×
359
360	const auto barcode = (barcode_it == barcode_end) ? 0 : *barcode_it++;	×
361	const auto alignment =	×
362	aligners.at(barcode).align_single_end(read, m_config.shift);	×
363
364	if (m_config.is_good_alignment(alignment)) {	×
365	const auto length = read.length();	×
366	alignment.truncate_single_end(read);	×
367
368	stats->adapter_trimmed_reads.inc(alignment.adapter_id);	×
369	stats->adapter_trimmed_bases.inc(alignment.adapter_id,	×
370	length - read.length());	×
371	}
372
373	// Add (optional) user specified prefixes to read names
374	read.add_prefix_to_name(m_config.prefix_read_1);	×
375
376	// Trim fixed number of bases from 5' and/or 3' termini
377	post_trim_read_termini(m_config, *stats, read, read_type::mate_1);	×
378	// Trim poly-X tails for zero or more X
379	post_trim_poly_x_tail(m_config, *stats, read);	×
380	// Sliding window trimming or single-base trimming of low quality bases
381	post_trim_read_by_quality(m_config, *stats, read);	×
382
383	stats->total_trimmed.inc_reads(in_length != read.length());	×
384	stats->total_trimmed.inc_bases(in_length - read.length());	×
385
386	if (is_acceptable_read(m_config, *stats, read)) {	×
387	stats->read_1->process(read);	×
388	chunks.add(read, read_type::mate_1, fastq_flags::se, barcode);	×
389	} else {
390	stats->discarded->process(read);	×
391	chunks.add(read, read_type::discarded, fastq_flags::se_fail, barcode);	×
392	}
393	}
394
395	m_stats.release(stats);	×
396
397	return chunks.finalize(chunk->eof);	×
398	}
399
400	////////////////////////////////////////////////////////////////////////////////
401	// Implementations for `pe_reads_processor`
402
403	void
404	add_pe_statistics(const trimming_statistics& stats,	×
405	const fastq& read,
406	read_type type)
407	{
408	switch (type) {	×
409	case read_type::mate_1:	×
410	stats.read_1->process(read);	×
411	break;	×
412	case read_type::mate_2:	×
413	stats.read_2->process(read);	×
414	break;	×
415
416	case read_type::merged:	×
417	stats.merged->process(read);	×
418	break;	×
419
420	case read_type::singleton:	×
421	stats.singleton->process(read);	×
422	break;	×
423
424	case read_type::discarded:	×
425	stats.discarded->process(read);	×
426	break;	×
427
428	case read_type::max:	×
429	AR_FAIL("unhandled read type");	×
430
431	default:	×
432	AR_FAIL("invalid read type");	×
433	}
434	}
435
436	pe_reads_processor::pe_reads_processor(const userconfig& config,	×
437	const sample_output_files& output,
438	const size_t nth,
439	trim_stats_ptr sink)	×
440	: reads_processor(config, output, nth, sink)	×
441	{
442	}
443
444	chunk_vec
445	pe_reads_processor::process(chunk_ptr chunk)	×
446	{
447	AR_REQUIRE(chunk);	×
448	processed_reads chunks{ m_output };	×
449	chunks.set_sample(m_config.samples.at(m_sample));	×
450	chunks.set_mate_separator(chunk->mate_separator);	×
451
452	if (chunk->first) {	×
453	chunks.write_headers(m_config.args);	×
454	}
455
456	sequence_merger merger;	×
457	merger.set_merge_strategy(m_config.merge);	×
458	merger.set_max_recalculated_score(m_config.merge_quality_max);	×
459
460	// A sequence aligner per barcode (pair)
461	std::vector<sequence_aligner> aligners;	×
462	for (const auto& it : m_config.samples.at(m_sample)) {	×
463	aligners.emplace_back(it.adapters, m_config.simd);	×
464	}
465
466	auto stats = m_stats.acquire();	×
467	stats->adapter_trimmed_reads.resize_up_to(m_config.samples.adapters().size());	×
468	stats->adapter_trimmed_bases.resize_up_to(m_config.samples.adapters().size());	×
469
470	AR_REQUIRE(!aligners.empty());	×
471	AR_REQUIRE(chunk->reads_1.size() == chunk->reads_2.size());	×
472	AR_REQUIRE(chunk->barcodes.empty() \|\|	×
473	chunk->barcodes.size() == chunk->reads_1.size());
474
475	auto barcode_it = chunk->barcodes.begin();	×
476	const auto barcode_end = chunk->barcodes.end();	×
477
478	auto it_1 = chunk->reads_1.begin();	×
479	auto it_2 = chunk->reads_2.begin();	×
480	while (it_1 != chunk->reads_1.end()) {	×
481	fastq& read_1 = *it_1++;	×
482	fastq& read_2 = *it_2++;	×
483
484	const size_t in_length_1 = read_1.length();	×
485	const size_t in_length_2 = read_2.length();	×
486
487	// Trim fixed number of bases from 5' and/or 3' termini
488	pre_trim_read_termini(m_config, *stats, read_1, read_type::mate_1);	×
489	pre_trim_read_termini(m_config, *stats, read_2, read_type::mate_2);	×
490
491	// Trim poly-X tails for zero or more X
492	pre_trim_poly_x_tail(m_config, *stats, read_1);	×
493	pre_trim_poly_x_tail(m_config, *stats, read_2);	×
494
495	// Reverse complement to match the orientation of read_1
496	read_2.reverse_complement();	×
497
498	const auto barcode = (barcode_it == barcode_end) ? 0 : *barcode_it++;	×
499	const auto alignment =	×
500	aligners.at(barcode).align_paired_end(read_1, read_2, m_config.shift);	×
501
502	if (m_config.is_good_alignment(alignment)) {	×
503	const size_t insert_size = alignment.insert_size(read_1, read_2);	×
504	const bool can_merge_alignment = m_config.can_merge_alignment(alignment);	×
505	if (can_merge_alignment) {	×
506	// Insert size calculated from untrimmed reads
507	stats->insert_sizes.resize_up_to(insert_size + 1);	×
508	stats->insert_sizes.inc(insert_size);	×
509	}
510
511	const size_t pre_trimmed_bp = read_1.length() + read_2.length();	×
512	const size_t n_adapters = alignment.truncate_paired_end(read_1, read_2);	×
513	const size_t post_trimmed_bp = read_1.length() + read_2.length();	×
514
515	stats->adapter_trimmed_reads.inc(alignment.adapter_id, n_adapters);	×
516	stats->adapter_trimmed_bases.inc(alignment.adapter_id,	×
517	pre_trimmed_bp - post_trimmed_bp);	×
518
519	if (m_config.merge != merge_strategy::none && can_merge_alignment) {	×
520	// Track if one or both source reads are trimmed post merging
521	merged_reads mstats(read_1, read_2, insert_size);	×
522
523	// Merge read_2 into read_1
524	merger.merge(alignment, read_1, read_2);	×
525
526	// Add (optional) user specified prefix to read names
527	read_1.add_prefix_to_name(m_config.prefix_merged);	×
528
529	// Trim fixed number of bases from 5' and/or 3' termini
530	post_trim_read_termini(	×
531	m_config, *stats, read_1, read_type::merged, &mstats);	×
532
533	if (!m_config.preserve5p) {	×
534	// A merged read essentially consists of two 5p termini, so neither
535	// end should be trimmed if --preserve5p is used.
536	post_trim_read_by_quality(m_config, *stats, read_1, &mstats);	×
537	}
538
539	// Track total amount of bases trimmed/lost
540	stats->total_trimmed.inc_reads(2);	×
541	stats->total_trimmed.inc_bases(in_length_1 + in_length_2 -	×
542	read_1.length());	×
543
544	if (is_acceptable_read(m_config, *stats, read_1, 2)) {	×
545	stats->merged->process(read_1, 2);	×
546	chunks.add(read_1, read_type::merged, fastq_flags::se, barcode);	×
547	} else {
548	stats->discarded->process(read_1, 2);	×
549	chunks.add(	×
550	read_1, read_type::discarded, fastq_flags::se_fail, barcode);
551	}
552
553	continue;	×
554	}
555	}
556
557	// Reads were not aligned or merging is not enabled
558	// Undo reverse complementation (post truncation of adapters)
559	read_2.reverse_complement();	×
560
561	// Add (optional) user specified prefixes to read names
562	read_1.add_prefix_to_name(m_config.prefix_read_1);	×
563	read_2.add_prefix_to_name(m_config.prefix_read_2);	×
564
565	// Trim fixed number of bases from 5' and/or 3' termini
566	post_trim_read_termini(m_config, *stats, read_1, read_type::mate_1);	×
567	post_trim_read_termini(m_config, *stats, read_2, read_type::mate_2);	×
568
569	// Trim poly-X tails for zero or more X
570	post_trim_poly_x_tail(m_config, *stats, read_1);	×
571	post_trim_poly_x_tail(m_config, *stats, read_2);	×
572
573	// Sliding window trimming or single-base trimming of low quality bases
574	post_trim_read_by_quality(m_config, *stats, read_1);	×
575	post_trim_read_by_quality(m_config, *stats, read_2);	×
576
577	// Are the reads good enough? Not too many Ns?
578	const bool is_ok_1 = is_acceptable_read(m_config, *stats, read_1);	×
579	const bool is_ok_2 = is_acceptable_read(m_config, *stats, read_2);	×
580
581	read_type type_1;	×
582	read_type type_2;	×
583	if (is_ok_1 && is_ok_2) {	×
584	type_1 = read_type::mate_1;
585	type_2 = read_type::mate_2;
586	} else if (is_ok_1) {	×
587	type_1 = read_type::singleton;
588	type_2 = read_type::discarded;
589	} else if (is_ok_2) {	×
590	type_1 = read_type::discarded;
591	type_2 = read_type::singleton;
592	} else {
593	type_1 = read_type::discarded;	×
594	type_2 = read_type::discarded;	×
595	}
596
597	add_pe_statistics(*stats, read_1, type_1);	×
598	add_pe_statistics(*stats, read_2, type_2);	×
599
600	// Total number of reads/bases trimmed
601	stats->total_trimmed.inc_reads((in_length_1 != read_1.length()) +	×
602	(in_length_2 != read_2.length()));	×
603	stats->total_trimmed.inc_bases((in_length_1 - read_1.length()) +	×
604	(in_length_2 - read_2.length()));	×
605
606	const auto flags_1 = is_ok_1 ? fastq_flags::pe_1 : fastq_flags::pe_1_fail;	×
607	const auto flags_2 = is_ok_2 ? fastq_flags::pe_2 : fastq_flags::pe_2_fail;	×
608
609	// Queue reads last, since this result in modifications to lengths
610	chunks.add(read_1, type_1, flags_1, barcode);	×
611	chunks.add(read_2, type_2, flags_2, barcode);	×
612	}
613
614	// Track amount of overlapping bases "lost" due to read merging
615	stats->reads_merged.inc_reads(merger.reads_merged());	×
616	stats->reads_merged.inc_bases(merger.bases_merged());	×
617
618	m_stats.release(stats);	×
619
620	return chunks.finalize(chunk->eof);	×
621	}
622
623	} // namespace adapterremoval	×

MikkelSchubert / adapterremoval / #45

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous