#102

Committed 18 Apr 2025 01:14PM UTC coverage: 67.126% (-0.06%) from 67.186%

Build # #102

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

support normalizing (merged) reverse inserts (#128)

This implements part of #68

Run Details

0 of 16 new or added lines in 3 files covered. (0.0%)

3 existing lines in 1 file now uncovered.

9697 of 14446 relevant lines covered (67.13%)

3061.02 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, ...
#include "commontypes.hpp"   // for read_file, trimming_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, ...
#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_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"
5	#include "alignment.hpp" // for alignment_info, sequence_merger, ...
6	#include "commontypes.hpp" // for read_file, trimming_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, ...
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
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;	×
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) {	×
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
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
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
171	case read_file::singleton:	×
172	case read_file::discarded:	×
173	AR_FAIL("unsupported read type in post_trim_read_termini");	×
174
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
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
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);	×
372	chunks.add(read, read_type::se, barcode);	×
373	} else {
374	stats->discarded->process(read);	×
375	chunks.add(read, 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	namespace {
388
389	void
390	add_pe_statistics(const trimming_statistics& stats,	×
391	const fastq& read,
392	read_file type)
393	{
394	switch (type) {	×
395	case read_file::mate_1:	×
396	stats.read_1->process(read);	×
397	break;	×
398	case read_file::mate_2:	×
399	stats.read_2->process(read);	×
400	break;	×
401
402	case read_file::singleton:	×
403	stats.singleton->process(read);	×
404	break;	×
405
406	case read_file::discarded:	×
407	stats.discarded->process(read);	×
408	break;	×
409
410	case read_file::max:	×
411	case read_file::merged:	×
412	AR_FAIL("unhandled read type");	×
413
414	default:	×
415	AR_FAIL("invalid read type");	×
416	}
417	}
418
419	} // namespace
420
421	pe_reads_processor::pe_reads_processor(const userconfig& config,	×
422	const sample_output_files& output,
423	const size_t nth,
424	trim_stats_ptr sink)	×
425	: reads_processor(config, output, nth, sink)	×
426	{
427	}
428
429	chunk_vec
430	pe_reads_processor::process(chunk_ptr chunk)	×
431	{
432	AR_REQUIRE(chunk);	×
433	processed_reads chunks{ m_output };	×
434	chunks.set_sample(m_config.samples.at(m_sample));	×
435	chunks.set_mate_separator(chunk->mate_separator);	×
436
437	if (chunk->first) {	×
438	chunks.write_headers(m_config.args);	×
439	}
440
441	sequence_merger merger;	×
442	merger.set_merge_strategy(m_config.merge);	×
443	merger.set_max_recalculated_score(m_config.merge_quality_max);	×
444
445	// A sequence aligner per barcode (pair)
446	std::vector<sequence_aligner> aligners;	×
NEW 447	const auto& sample = m_config.samples.at(m_sample);	×
NEW 448	for (const auto& it : sample) {	×
UNCOV 449	aligners.emplace_back(it.adapters, m_config.simd);	×
450	}
451
452	auto stats = m_stats.acquire();	×
453	stats->adapter_trimmed_reads.resize_up_to(m_config.samples.adapters().size());	×
454	stats->adapter_trimmed_bases.resize_up_to(m_config.samples.adapters().size());	×
455
456	AR_REQUIRE(!aligners.empty());	×
457	AR_REQUIRE(chunk->reads_1.size() == chunk->reads_2.size());	×
458	AR_REQUIRE(chunk->barcodes.empty() \|\|	×
459	chunk->barcodes.size() == chunk->reads_1.size());
460
461	auto barcode_it = chunk->barcodes.begin();	×
462	const auto barcode_end = chunk->barcodes.end();	×
463
464	auto it_1 = chunk->reads_1.begin();	×
465	auto it_2 = chunk->reads_2.begin();	×
466	while (it_1 != chunk->reads_1.end()) {	×
467	fastq& read_1 = *it_1++;	×
468	fastq& read_2 = *it_2++;	×
469
470	const size_t in_length_1 = read_1.length();	×
471	const size_t in_length_2 = read_2.length();	×
472
473	// Trim fixed number of bases from 5' and/or 3' termini
474	pre_trim_read_termini(m_config, *stats, read_1, read_file::mate_1);	×
475	pre_trim_read_termini(m_config, *stats, read_2, read_file::mate_2);	×
476
477	// Trim poly-X tails for zero or more X
478	pre_trim_poly_x_tail(m_config, *stats, read_1);	×
479	pre_trim_poly_x_tail(m_config, *stats, read_2);	×
480
481	// Reverse complement to match the orientation of read_1
482	read_2.reverse_complement();	×
483
484	const auto barcode = (barcode_it == barcode_end) ? 0 : *barcode_it++;	×
485	const auto alignment =	×
486	aligners.at(barcode).align_paired_end(read_1, read_2, m_config.shift);	×
487
488	if (m_config.is_good_alignment(alignment)) {	×
489	const size_t insert_size = alignment.insert_size(read_1, read_2);	×
490	const bool can_merge_alignment = m_config.can_merge_alignment(alignment);	×
491	if (can_merge_alignment) {	×
492	// Insert size calculated from untrimmed reads
493	stats->insert_sizes.resize_up_to(insert_size + 1);	×
494	stats->insert_sizes.inc(insert_size);	×
495	}
496
497	const size_t pre_trimmed_bp = read_1.length() + read_2.length();	×
498	const size_t n_adapters = alignment.truncate_paired_end(read_1, read_2);	×
499	const size_t post_trimmed_bp = read_1.length() + read_2.length();	×
500
501	stats->adapter_trimmed_reads.inc(alignment.adapter_id, n_adapters);	×
502	stats->adapter_trimmed_bases.inc(alignment.adapter_id,	×
503	pre_trimmed_bp - post_trimmed_bp);	×
504
505	if (m_config.merge != merge_strategy::none && can_merge_alignment) {	×
506	// Track if one or both source reads are trimmed post merging
507	merged_reads mstats(read_1, read_2, insert_size);	×
508
509	// Merge read_2 into read_1
510	merger.merge(alignment, read_1, read_2);	×
511
512	// Add (optional) user specified prefix to read names
513	read_1.add_prefix_to_name(m_config.prefix_merged);	×
514
515	// Trim fixed number of bases from 5' and/or 3' termini
516	post_trim_read_termini(m_config,	×
517	*stats,	×
518	read_1,
519	read_file::merged,
520	&mstats);
521
522	if (!m_config.preserve5p) {	×
523	// A merged read essentially consists of two 5p termini, so neither
524	// end should be trimmed if --preserve5p is used.
525	post_trim_read_by_quality(m_config, *stats, read_1, &mstats);	×
526	}
527
528	// Track total amount of bases trimmed/lost
529	stats->total_trimmed.inc_reads(2);	×
530	stats->total_trimmed.inc_bases(in_length_1 + in_length_2 -	×
531	read_1.length());	×
532
NEW 533	auto meta = read_meta(read_type::merged).barcode(barcode);	×
NEW 534	if (m_config.normalize_orientation &&	×
NEW 535	sample.at(barcode).orientation == barcode_orientation::reverse) {	×
NEW 536	read_1.reverse_complement();	×
537	}
538
539	if (is_acceptable_read(m_config, *stats, read_1, 2)) {	×
540	stats->merged->process(read_1, 2);	×
541	} else {
UNCOV 542	stats->discarded->process(read_1, 2);	×
NEW 543	meta.type(read_type::merged_fail);	×
544	}
545
NEW 546	chunks.add(read_1, meta);	×
547
UNCOV 548	continue;	×
549	}
550	}
551
552	// Reads were not aligned or merging is not enabled
553	// Undo reverse complementation (post truncation of adapters)
554	read_2.reverse_complement();	×
555
556	// Add (optional) user specified prefixes to read names
557	read_1.add_prefix_to_name(m_config.prefix_read_1);	×
558	read_2.add_prefix_to_name(m_config.prefix_read_2);	×
559
560	// Trim fixed number of bases from 5' and/or 3' termini
561	post_trim_read_termini(m_config, *stats, read_1, read_file::mate_1);	×
562	post_trim_read_termini(m_config, *stats, read_2, read_file::mate_2);	×
563
564	// Trim poly-X tails for zero or more X
565	post_trim_poly_x_tail(m_config, *stats, read_1);	×
566	post_trim_poly_x_tail(m_config, *stats, read_2);	×
567
568	// Sliding window trimming or single-base trimming of low quality bases
569	post_trim_read_by_quality(m_config, *stats, read_1);	×
570	post_trim_read_by_quality(m_config, *stats, read_2);	×
571
572	// Are the reads good enough? Not too many Ns?
573	const bool is_ok_1 = is_acceptable_read(m_config, *stats, read_1);	×
574	const bool is_ok_2 = is_acceptable_read(m_config, *stats, read_2);	×
575
576	read_meta meta_1{ read_type::pe_1 };	×
577	read_meta meta_2{ read_type::pe_2 };	×
578	if (!is_ok_1 \|\| !is_ok_2) {	×
579	if (is_ok_1) {	×
580	meta_1 = read_meta{ read_type::singleton_1 };	×
581	meta_2 = read_meta{ read_type::pe_2_fail };	×
582	} else if (is_ok_2) {	×
583	meta_1 = read_meta{ read_type::pe_1_fail };	×
584	meta_2 = read_meta{ read_type::singleton_2 };	×
585	} else {
586	meta_1 = read_meta{ read_type::pe_1_fail };	×
587	meta_2 = read_meta{ read_type::pe_2_fail };	×
588	}
589	}
590
591	add_pe_statistics(*stats, read_1, meta_1.get_file());	×
592	add_pe_statistics(*stats, read_2, meta_2.get_file());	×
593
594	// Total number of reads/bases trimmed
595	stats->total_trimmed.inc_reads((in_length_1 != read_1.length()) +	×
596	(in_length_2 != read_2.length()));	×
597	stats->total_trimmed.inc_bases((in_length_1 - read_1.length()) +	×
598	(in_length_2 - read_2.length()));	×
599
600	// Queue reads last, since this result in modifications to lengths
601	chunks.add(read_1, meta_1.barcode(barcode));	×
602	chunks.add(read_2, meta_2.barcode(barcode));	×
603	}
604
605	// Track amount of overlapping bases "lost" due to read merging
606	stats->reads_merged.inc_reads(merger.reads_merged());	×
607	stats->reads_merged.inc_bases(merger.bases_merged());	×
608
609	m_stats.release(stats);	×
610
611	return chunks.finalize(chunk->eof);	×
612	}
613
614	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #102

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