MikkelSchubert / adapterremoval / #95

Committed 13 Apr 2025 06:28PM UTC coverage: 30.749%. Remained the same

Build # #95

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

rename enums read_type and fastq_flags (#121)

To clarifies usage of two commonly used enums, `read_type` was renamed
to `read_file`, since it is used to specify input and output file
types, while `fastq_flags` was renamed to `read_type` since it is used
to differentiate between different pre/post-processed read types

Run Details

0 of 66 new or added lines in 6 files covered. (0.0%)

4 existing lines in 1 file now uncovered.

3157 of 10267 relevant lines covered (30.75%)

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

        if (is_acceptable_read(m_config, *stats, read_1, 2)) {
          stats->merged->process(read_1, 2);
          chunks.add(read_1, read_file::merged, read_type::se, barcode);
        } else {
          stats->discarded->process(read_1, 2);
          chunks.add(read_1, read_file::discarded, read_type::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_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_file type_1;
    read_file type_2;
    if (is_ok_1 && is_ok_2) {
      type_1 = read_file::mate_1;
      type_2 = read_file::mate_2;
    } else if (is_ok_1) {
      type_1 = read_file::singleton;
      type_2 = read_file::discarded;
    } else if (is_ok_2) {
      type_1 = read_file::discarded;
      type_2 = read_file::singleton;
    } else {
      type_1 = read_file::discarded;
      type_2 = read_file::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 ? read_type::pe_1 : read_type::pe_1_fail;
    const auto flags_2 = is_ok_2 ? read_type::pe_2 : read_type::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_file, 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 read_type
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_file type)
122	{
123	size_t trim_5p = 0;	×
124	size_t trim_3p = 0;	×
125
NEW 126	if (type == read_file::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;	×
NEW 131	} else if (type == read_file::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_file type,
149	merged_reads* mstats = nullptr)
150	{
151	size_t trim_5p = 0;	×
152	size_t trim_3p = 0;	×
153
154	switch (type) {	×
NEW 155	case read_file::mate_1:	×
156	trim_5p = config.post_trim_fixed_5p.first;	×
157	trim_3p = config.post_trim_fixed_3p.first;	×
158	break;	×
159
NEW 160	case read_file::mate_2:	×
161	trim_5p = config.post_trim_fixed_5p.second;	×
162	trim_3p = config.post_trim_fixed_3p.second;	×
163	break;	×
164
NEW 165	case read_file::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
NEW 171	case read_file::singleton:	×
NEW 172	case read_file::discarded:	×
UNCOV 173	AR_FAIL("unsupported read type in post_trim_read_termini");	×
174
NEW 175	case read_file::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
NEW 340	pre_trim_read_termini(m_config, *stats, read, read_file::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
NEW 361	post_trim_read_termini(m_config, *stats, read, read_file::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);	×
NEW 372	chunks.add(read, read_file::mate_1, read_type::se, barcode);	×
373	} else {
374	stats->discarded->process(read);	×
NEW 375	chunks.add(read, read_file::discarded, read_type::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_file type)
391	{
392	switch (type) {	×
NEW 393	case read_file::mate_1:	×
394	stats.read_1->process(read);	×
395	break;	×
NEW 396	case read_file::mate_2:	×
397	stats.read_2->process(read);	×
398	break;	×
399
NEW 400	case read_file::merged:	×
401	stats.merged->process(read);	×
402	break;	×
403
NEW 404	case read_file::singleton:	×
405	stats.singleton->process(read);	×
406	break;	×
407
NEW 408	case read_file::discarded:	×
409	stats.discarded->process(read);	×
410	break;	×
411
NEW 412	case read_file::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
NEW 472	pre_trim_read_termini(m_config, *stats, read_1, read_file::mate_1);	×
NEW 473	pre_trim_read_termini(m_config, *stats, read_2, read_file::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_file::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);	×
NEW 533	chunks.add(read_1, read_file::merged, read_type::se, barcode);	×
534	} else {
535	stats->discarded->process(read_1, 2);	×
NEW 536	chunks.add(read_1, read_file::discarded, read_type::se_fail, barcode);	×
537	}
538
539	continue;	×
540	}
541	}
542
543	// Reads were not aligned or merging is not enabled
544	// Undo reverse complementation (post truncation of adapters)
545	read_2.reverse_complement();	×
546
547	// Add (optional) user specified prefixes to read names
548	read_1.add_prefix_to_name(m_config.prefix_read_1);	×
549	read_2.add_prefix_to_name(m_config.prefix_read_2);	×
550
551	// Trim fixed number of bases from 5' and/or 3' termini
NEW 552	post_trim_read_termini(m_config, *stats, read_1, read_file::mate_1);	×
NEW 553	post_trim_read_termini(m_config, *stats, read_2, read_file::mate_2);	×
554
555	// Trim poly-X tails for zero or more X
556	post_trim_poly_x_tail(m_config, *stats, read_1);	×
557	post_trim_poly_x_tail(m_config, *stats, read_2);	×
558
559	// Sliding window trimming or single-base trimming of low quality bases
560	post_trim_read_by_quality(m_config, *stats, read_1);	×
561	post_trim_read_by_quality(m_config, *stats, read_2);	×
562
563	// Are the reads good enough? Not too many Ns?
564	const bool is_ok_1 = is_acceptable_read(m_config, *stats, read_1);	×
565	const bool is_ok_2 = is_acceptable_read(m_config, *stats, read_2);	×
566
NEW 567	read_file type_1;	×
NEW 568	read_file type_2;	×
UNCOV 569	if (is_ok_1 && is_ok_2) {	×
570	type_1 = read_file::mate_1;
571	type_2 = read_file::mate_2;
UNCOV 572	} else if (is_ok_1) {	×
573	type_1 = read_file::singleton;
574	type_2 = read_file::discarded;
UNCOV 575	} else if (is_ok_2) {	×
576	type_1 = read_file::discarded;
577	type_2 = read_file::singleton;
578	} else {
NEW 579	type_1 = read_file::discarded;	×
NEW 580	type_2 = read_file::discarded;	×
581	}
582
583	add_pe_statistics(*stats, read_1, type_1);	×
584	add_pe_statistics(*stats, read_2, type_2);	×
585
586	// Total number of reads/bases trimmed
587	stats->total_trimmed.inc_reads((in_length_1 != read_1.length()) +	×
588	(in_length_2 != read_2.length()));	×
589	stats->total_trimmed.inc_bases((in_length_1 - read_1.length()) +	×
590	(in_length_2 - read_2.length()));	×
591
NEW 592	const auto flags_1 = is_ok_1 ? read_type::pe_1 : read_type::pe_1_fail;	×
NEW 593	const auto flags_2 = is_ok_2 ? read_type::pe_2 : read_type::pe_2_fail;	×
594
595	// Queue reads last, since this result in modifications to lengths
596	chunks.add(read_1, type_1, flags_1, barcode);	×
597	chunks.add(read_2, type_2, flags_2, barcode);	×
598	}
599
600	// Track amount of overlapping bases "lost" due to read merging
601	stats->reads_merged.inc_reads(merger.reads_merged());	×
602	stats->reads_merged.inc_bases(merger.bases_merged());	×
603
604	m_stats.release(stats);	×
605
606	return chunks.finalize(chunk->eof);	×
607	}
608
609	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #95

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