MikkelSchubert / adapterremoval / #78

Committed 03 Apr 2025 03:38PM UTC coverage: 27.703% (-0.08%) from 27.778%

Build # #78

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

improve ccache hit rates (#101)

There is generally a high amount of churn in ccache files, due to the (currently) high level of interdependence between files, and cached files are therefore not expected have be relevant for a long time. The maximum size of the cache is therefore reduced to 100 MB, compared to the default maximum of 500 MB, in order to reduce the amount of time spent on (re)storing the cache.

Additionally, some tweaks were made to the build process, to ensure that building was separated from installing and testing, and `ubuntu-22.04` was replaced with `ubuntu-24.04` in order to (hopefully) fix an issue were unit test compilation was not ccache'd

Run Details

2701 of 9750 relevant lines covered (27.7%)

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

MikkelSchubert / adapterremoval / #78

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