#117

Committed 25 May 2025 03:01PM UTC coverage: 66.932% (-0.07%) from 67.006%

Build # #117

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

iwyu and reduce build-time inter-dependencies (#144)

Run Details

26 of 145 new or added lines in 20 files covered. (17.93%)

89 existing lines in 5 files now uncovered.

9738 of 14549 relevant lines covered (66.93%)

3041.19 hits per line

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

0.0

/src/trimming.cpp

// SPDX-License-Identifier: GPL-3.0-or-later
// SPDX-FileCopyrightText: 2011 Stinus Lindgreen <stinus@binf.ku.dk>
// SPDX-FileCopyrightText: 2014 Mikkel Schubert <mikkelsch@gmail.com>
#include "trimming.hpp"      // declarations
#include "alignment.hpp"     // for alignment_info, sequence_merger, ...
#include "commontypes.hpp"   // for read_type, trimming_strategy, ...
#include "counts.hpp"        // for counts, indexed_count
#include "debug.hpp"         // for AR_FAIL, AR_REQUIRE
#include "fastq.hpp"         // for fastq
#include "fastq_enc.hpp"     // for ACGT
#include "output.hpp"        // for sample_output_files, processed_reads
#include "scheduler.hpp"     // for analytical_step, processing_order
#include "sequence_sets.hpp" // for adapter_set
#include "serializer.hpp"    // for read_type
#include "statistics.hpp"    // for trimming_statistics, reads_and_bases, ...
#include "userconfig.hpp"    // for userconfig
#include <memory>            // for unique_ptr, __shared_ptr_access, make_unique
#include <string>            // for string
#include <utility>           // for pair, move
#include <vector>            // for vector

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_file type)
{
  size_t trim_5p = 0;
  size_t trim_3p = 0;

  if (type == read_file::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_file::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_file type,
                       merged_reads* mstats = nullptr)
{
  size_t trim_5p = 0;
  size_t trim_3p = 0;

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

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

    case read_file::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_file::singleton:
    case read_file::discarded:
      AR_FAIL("unsupported read type in post_trim_read_termini");

    case read_file::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_file::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_file::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::se, barcode);
    } else {
      stats->discarded->process(read);
      chunks.add(read, read_type::se_fail, barcode);
    }
  }

  m_stats.release(stats);

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

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

namespace {

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

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

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

    case read_file::max:
    case read_file::merged:
      AR_FAIL("unhandled read type");

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

} // namespace

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;
  const auto& sample = m_config.samples->at(m_sample);
  for (const auto& it : 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_file::mate_1);
    pre_trim_read_termini(m_config, *stats, read_2, read_file::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_file::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());

        auto meta = read_meta(read_type::merged).barcode(barcode);
        if (m_config.normalize_orientation &&
            sample.at(barcode).orientation == barcode_orientation::reverse) {
          read_1.reverse_complement();
        }

        if (is_acceptable_read(m_config, *stats, read_1, 2)) {
          stats->merged->process(read_1, 2);
        } else {
          stats->discarded->process(read_1, 2);
          meta.type(read_type::merged_fail);
        }

        chunks.add(read_1, meta);

        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_file::mate_1);
    post_trim_read_termini(m_config, *stats, read_2, read_file::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_meta meta_1{ read_type::pe_1 };
    read_meta meta_2{ read_type::pe_2 };
    if (!is_ok_1 || !is_ok_2) {
      if (is_ok_1) {
        meta_1 = read_meta{ read_type::singleton_1 };
        meta_2 = read_meta{ read_type::pe_2_fail };
      } else if (is_ok_2) {
        meta_1 = read_meta{ read_type::pe_1_fail };
        meta_2 = read_meta{ read_type::singleton_2 };
      } else {
        meta_1 = read_meta{ read_type::pe_1_fail };
        meta_2 = read_meta{ read_type::pe_2_fail };
      }
    }

    add_pe_statistics(*stats, read_1, meta_1.get_file());
    add_pe_statistics(*stats, read_2, meta_2.get_file());

    // 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()));

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

MikkelSchubert / adapterremoval / #117

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