MikkelSchubert / adapterremoval / #72

Committed 22 Mar 2025 08:35PM UTC coverage: 27.09% (+0.02%) from 27.075%

Build # #72

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

rework unit tests to check exception messages (#94)

This includes fixing a pointless test and a improper error messages on
lowercase degenerate FASTQ bases

Coverage Stats

4 of 5 new or added lines in 1 file covered. (80.0%)

2599 of 9594 relevant lines covered (27.09%)

4262.56 hits per line

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

75.8

/src/fastq_enc.cpp

/*************************************************************************\
 * AdapterRemoval - cleaning next-generation sequencing reads            *
 *                                                                       *
 * Copyright (C) 2011 by Stinus Lindgreen - stinus@binf.ku.dk            *
 * Copyright (C) 2014 by Mikkel Schubert - mikkelsch@gmail.com           *
 *                                                                       *
 * This program is free software: you can redistribute it and/or modify  *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * This program is distributed in the hope that it will be useful,       *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have received a copy of the GNU General Public License     *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>. *
\*************************************************************************/
#include "fastq_enc.hpp" // header
#include "debug.hpp"     // for AR_FAIL, AR_REQUIRE
#include "errors.hpp"    // for fastq_error
#include "strutils.hpp"  // for shell_escape
#include <algorithm>     // for min, max
#include <array>         // for array
#include <cmath>         // for log10, pow, round
#include <sstream>       // for ostringstream

namespace adapterremoval {

namespace {

//! Offset used by Phred+33 and SAM encodings
const int PHRED_33_OFFSET_MIN = '!';
//! The maximum ASCII value allowed for Phred+33 encoded quality scores. This is
//! a generous maximum, taking into account that instruments are assigning
//! higher and higher quality scores for Phred+33 data.
const int PHRED_33_OFFSET_MAX = 'N';
//! Minimum Phred score allowed
// const int PHRED_33_SCORE_MIN = 0;
//! Maximum Phred score allowed
const int PHRED_33_SCORE_MAX = PHRED_33_OFFSET_MAX - PHRED_33_OFFSET_MIN;

//! Offset used by Phred+64 encodings
const int PHRED_64_OFFSET_MIN = '@';
//! The maximum ASCII value allowed for encoded Phred scores
const int PHRED_64_OFFSET_MAX = '~';
//! Minimum Phred+64 score allowed
const int PHRED_64_SCORE_MIN = 0;
//! Maximum Phred+64 score allowed
// const int PHRED_64_SCORE_MAX = PHRED_64_OFFSET_MAX - PHRED_64_OFFSET_MIN;

//! Offset used by Solexa encoding quality scores
const int SOLEXA_OFFSET_MIN = '@';
//! The maximum ASCII value allowed for encoded Solexa scores
const int SOLEXA_OFFSET_MAX = 'h';
//! Minimum Phred encoded score allowed
const int SOLEXA_SCORE_MIN = -5;
//! Maximum Phred encoded score allowed
const int SOLEXA_SCORE_MAX = SOLEXA_OFFSET_MAX - SOLEXA_OFFSET_MIN;

///////////////////////////////////////////////////////////////////////////////
// Pre-calculation of Solexa <-> Phred conversions

std::array<char, 256>
calc_solexa_to_phred33()
{
  std::array<char, 256> scores = { 0 };

  for (int i = SOLEXA_SCORE_MIN; i <= SOLEXA_SCORE_MAX; ++i) {
    const double score = round(10.0 * log10(1.0 + pow(10, (i / 10.0))));
    const int transformed =
      std::max<int>(PHRED_SCORE_MIN, std::min<int>(PHRED_SCORE_MAX, score));
    scores.at(i + SOLEXA_OFFSET_MIN) = PHRED_OFFSET_MIN + transformed;
  }

  return scores;
}

const auto g_solexa_to_phred33 = calc_solexa_to_phred33();

///////////////////////////////////////////////////////////////////////////////

/**
 * Uppercase letters in the range a-z, but mangles other characters. Resulting
 * values outside of A-Z should be compared against/reported, since they may not
 * reflect the original value.
 */
constexpr char
mangle_to_upper(const char c)
{
  return c & 0xDF;
}

/** Returns the quality score in a form that is human readable */
std::string
escape_raw_score(char raw)
{
  return shell_escape(std::string(1, raw));
}

[[noreturn]] void
throw_invalid_phred_33(const char raw)
{
  const int score = raw - PHRED_33_OFFSET_MIN;
  const int alt_score = raw - PHRED_64_OFFSET_MIN;
  const int max_score = PHRED_33_SCORE_MAX;
  const auto esc = escape_raw_score(raw);
  const auto esc_max = escape_raw_score(PHRED_33_OFFSET_MAX);

  std::ostringstream ss;
  ss << "Found Phred+33 encoded quality score of " << score << " (encoded as "
     << esc << "), which is greater than the expected maximum score of "
     << max_score << " (encoded as " << esc_max << "). This suggests that "
     << "input may be Phred+64 encoded.\n\n"

     << "If the quality scores are actually Phred+64 encoded, which would "
     << "mean that the Phred quality score is " << alt_score
     << ", then use the '--quality-format 64' command-line option.\n\n"

     << "If the quality scores are Phred+33 encoded, but have higher than "
        "expected scores, then use '--quality-format sam' to permit the full "
        "range of quality scores.\n\n"

     << "See the documentation for more information.";

  throw fastq_error(ss.str());
}

[[noreturn]] void
throw_invalid_phred_64(const char raw)
{
  AR_REQUIRE(raw < SOLEXA_OFFSET_MIN,
             "invalid_phred called on valid PHRED score");

  const int score = raw - PHRED_64_OFFSET_MIN;
  const int alt_score = raw - PHRED_33_OFFSET_MIN;
  const int min_score = PHRED_64_SCORE_MIN;
  const auto esc = escape_raw_score(raw);
  const auto esc_min = escape_raw_score(PHRED_64_OFFSET_MIN);

  std::ostringstream ss;
  ss << "Found Phred+64 encoded quality score of " << score << " (encoded as "
     << esc << "), which is less than the expected minimum score of "
     << min_score << " (encoded as " << esc_min << "). This suggests that "
     << "input may be Phred+33 encoded.\n\n"

     << "If the quality scores are actually Phred+33 encoded, which would "
     << "mean that the Phred quality score is " << alt_score
     << ", then omit the '--quality-format' command-line option or use "
     << "'--quality-format 33' to explicitly set the format.\n\n";

  if (score >= SOLEXA_SCORE_MIN) {
    ss << "The quality score could also be the older Solexa format, which has "
       << "a minimum score of -5, but data of this type is rare. If it is "
       << "actually Solexa encoded FASTQ data, then use the '--quality-format "
       << "solexa' command-line option.\n\n";
  }

  ss << "See the documentation for more information.";

  throw fastq_error(ss.str());
}

[[noreturn]] void
throw_invalid_solexa(const char raw)
{
  const int score = raw - SOLEXA_OFFSET_MIN;
  const int alt_score = raw - PHRED_33_OFFSET_MIN;
  const int min_score = SOLEXA_SCORE_MIN;
  const auto esc = escape_raw_score(raw);
  const auto esc_min = escape_raw_score(SOLEXA_OFFSET_MIN + SOLEXA_SCORE_MIN);
  const auto esc_max = escape_raw_score(SOLEXA_OFFSET_MAX);

  std::ostringstream ss;
  if (score < SOLEXA_SCORE_MIN) {
    ss << "Found Solexa encoded quality score of " << score << " (encoded as "
       << esc << "), which is less than the expected minimum score of "
       << min_score << " (encoded as " << esc_min << "). This suggests that "
       << "input may be Phred+33 encoded.\n\n"

       << "If the quality scores are actually Phred+33 encoded, which would "
       << "mean that the Phred quality score is " << alt_score
       << ", then omit the '--quality-format' command-line option or use "
       << "'--quality-format 33' to explicitly set the format.\n\n"

       << "See the documentation for more information.";
  } else if (raw > SOLEXA_OFFSET_MAX) {
    ss << "Found Solexa encoded quality score of " << score << " (encoded as "
       << esc << "), which is greater than the expected maximum score of "
       << SOLEXA_SCORE_MAX << " (encoded as " << esc_max << ").\n\n"

       << "See the documentation for more information.";
  } else {
    AR_FAIL("invalid_phred called on valid PHRED score");
  }

  throw fastq_error(ss.str());
}

[[noreturn]] void
throw_invalid_score(const quality_encoding encoding, const char raw_score)
{
  if (raw_score < PHRED_OFFSET_MIN || raw_score > PHRED_OFFSET_MAX) {
    std::ostringstream ss;

    ss << "Found raw FASTQ quality score of " << static_cast<int>(raw_score)
       << " (" << escape_raw_score(raw_score) << "). This is outside of the "
       << "range of valid ASCII encoded quality scores (" << PHRED_OFFSET_MIN
       << " to " << PHRED_OFFSET_MAX << "), meaning that the input file is "
       << "either corrupt or not in FASTQ format!";

    throw fastq_error(ss.str());
  }

  switch (encoding) {
    case quality_encoding::phred_33:
      throw_invalid_phred_33(raw_score);

    case quality_encoding::phred_64:
      throw_invalid_phred_64(raw_score);

    case quality_encoding::solexa:
      throw_invalid_solexa(raw_score);

    case quality_encoding::sam:
      AR_FAIL("This case should have been handled by initial check");

    default:
      AR_FAIL("Invalid quality encoding");
  }
}

[[noreturn]] void
throw_invalid_base(char c)
{
  std::ostringstream stream;

  switch (mangle_to_upper(c)) {
    case 'B': // C / G / T
    case 'D': // A / G / T
    case 'H': // A / C / T
    case 'K': // G / T
    case 'M': // A / C
    case 'R': // A / G
    case 'S': // C / G
    case 'V': // A / C / G
    case 'W': // A / T
    case 'Y': // C / T
      stream << "found degenerate base '" << c << "' in FASTQ sequence, but "
             << "only bases A, C, G, T and N are supported. Use the option "
                "--mask-degenerate-bases to convert degenerate bases to N";
      break;

    case 'U': // Uracils
      stream << "found uracil (U) in FASTQ sequence, but only bases A, C, G, "
                "T and N are supported. Use the option --convert-uracils to "
                "convert uracils (U) to thymine (T)";
      break;

    default:
      stream << "invalid character " << log_escape(std::string(1, c))
             << " found in FASTQ sequence";
      break;
  }

  throw fastq_error(stream.str());
}

void
process_nucleotides_strict(std::string& nucleotides)
{
  for (char& nuc : nucleotides) {
    // Fast ASCII letter uppercase
    const auto upper_case = mangle_to_upper(nuc);
    switch (upper_case) {
      case 'A':
      case 'C':
      case 'G':
      case 'T':
      case 'N':
        nuc = upper_case;
        break;

      default:
        throw_invalid_base(nuc);
    }
  }
}

void
process_nucleotides_lenient(std::string& nucleotides,
                            const bool convert_uracil,
                            const bool mask_degenerate)
{
  for (char& nuc : nucleotides) {
    // Fast ASCII letter uppercase
    const auto upper_case = mangle_to_upper(nuc);
    switch (upper_case) {
      case 'A':
      case 'C':
      case 'G':
      case 'T':
      case 'N':
        nuc = upper_case;
        break;

      // Uracils
      case 'U':
        if (convert_uracil) {
          nuc = 'T';
          break;
        } else {
          throw_invalid_base(nuc);
        }

      // IUPAC encoded degenerate bases
      case 'B': // C / G / T
      case 'D': // A / G / T
      case 'H': // A / C / T
      case 'K': // G / T
      case 'M': // A / C
      case 'R': // A / G
      case 'S': // C / G
      case 'V': // A / C / G
      case 'W': // A / T
      case 'Y': // C / T
        if (mask_degenerate) {
          nuc = 'N';
          break;
        } else {
          throw_invalid_base(nuc);
        }

      default:
        throw_invalid_base(nuc);
    }
  }
}

} // namespace

///////////////////////////////////////////////////////////////////////////////

fastq_encoding::fastq_encoding(quality_encoding encoding,
                               degenerate_encoding degenerate,
                               uracil_encoding uracils) noexcept
  : m_mask_degenerate()
  , m_convert_uracil()
  , m_encoding(encoding)
  , m_offset_min()
  , m_offset_max()
{
  switch (degenerate) {
    case degenerate_encoding::reject:
      m_mask_degenerate = false;
      break;
    case degenerate_encoding::mask:
      m_mask_degenerate = true;
      break;
    default:
      AR_FAIL("invalid degenerate encoding value");
  }

  switch (uracils) {
    case uracil_encoding::convert:
      m_convert_uracil = true;
      break;
    case uracil_encoding::reject:
      m_convert_uracil = false;
      break;
    default:
      AR_FAIL("invalid uracil encoding value");
  }

  switch (encoding) {
    case quality_encoding::phred_33:
      m_offset_min = PHRED_33_OFFSET_MIN;
      m_offset_max = PHRED_33_OFFSET_MAX;
      break;

    case quality_encoding::phred_64:
      m_offset_min = PHRED_64_OFFSET_MIN;
      m_offset_max = PHRED_64_OFFSET_MAX;
      break;

    case quality_encoding::solexa:
      m_offset_min = SOLEXA_OFFSET_MIN;
      m_offset_max = SOLEXA_OFFSET_MAX;
      break;

    case quality_encoding::sam:
      m_offset_min = PHRED_OFFSET_MIN;
      m_offset_max = PHRED_OFFSET_MAX;
      break;

    default:
      AR_FAIL("unknown quality encoding");
  }
}

void
fastq_encoding::process_nucleotides(std::string& sequence) const
{
  if (m_mask_degenerate || m_convert_uracil) {
    process_nucleotides_lenient(sequence, m_convert_uracil, m_mask_degenerate);
  } else {
    process_nucleotides_strict(sequence);
  }
}

void
fastq_encoding::process_qualities(std::string& qualities) const
{
  if (m_encoding == quality_encoding::solexa) {
    for (auto& quality : qualities) {
      const char current = quality;
      // TODO: Handle negative values, e.g. è
      if (!(quality = g_solexa_to_phred33.at(quality))) {
        throw_invalid_score(m_encoding, current);
      }
    }
  } else {
    for (auto& quality : qualities) {
      if (quality < m_offset_min || quality > m_offset_max) {
        throw_invalid_score(m_encoding, quality);
      }

      // Convert Phred+64 to Phred+33 if needed
      quality -= m_offset_min - PHRED_33_OFFSET_MIN;
    }
  }
}

} // namespace adapterremoval

1	/*************************************************************************\
2	* AdapterRemoval - cleaning next-generation sequencing reads *
3	* *
4	* Copyright (C) 2011 by Stinus Lindgreen - stinus@binf.ku.dk *
5	* Copyright (C) 2014 by Mikkel Schubert - mikkelsch@gmail.com *
6	* *
7	* This program is free software: you can redistribute it and/or modify *
8	* it under the terms of the GNU General Public License as published by *
9	* the Free Software Foundation, either version 3 of the License, or *
10	* (at your option) any later version. *
11	* *
12	* This program is distributed in the hope that it will be useful, *
13	* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15	* GNU General Public License for more details. *
16	* *
17	* You should have received a copy of the GNU General Public License *
18	* along with this program. If not, see <http://www.gnu.org/licenses/>. *
19	\*************************************************************************/
20	#include "fastq_enc.hpp" // header
21	#include "debug.hpp" // for AR_FAIL, AR_REQUIRE
22	#include "errors.hpp" // for fastq_error
23	#include "strutils.hpp" // for shell_escape
24	#include <algorithm> // for min, max
25	#include <array> // for array
26	#include <cmath> // for log10, pow, round
27	#include <sstream> // for ostringstream
28
29	namespace adapterremoval {
30
31	namespace {
32
33	//! Offset used by Phred+33 and SAM encodings
34	const int PHRED_33_OFFSET_MIN = '!';
35	//! The maximum ASCII value allowed for Phred+33 encoded quality scores. This is
36	//! a generous maximum, taking into account that instruments are assigning
37	//! higher and higher quality scores for Phred+33 data.
38	const int PHRED_33_OFFSET_MAX = 'N';
39	//! Minimum Phred score allowed
40	// const int PHRED_33_SCORE_MIN = 0;
41	//! Maximum Phred score allowed
42	const int PHRED_33_SCORE_MAX = PHRED_33_OFFSET_MAX - PHRED_33_OFFSET_MIN;
43
44	//! Offset used by Phred+64 encodings
45	const int PHRED_64_OFFSET_MIN = '@';
46	//! The maximum ASCII value allowed for encoded Phred scores
47	const int PHRED_64_OFFSET_MAX = '~';
48	//! Minimum Phred+64 score allowed
49	const int PHRED_64_SCORE_MIN = 0;
50	//! Maximum Phred+64 score allowed
51	// const int PHRED_64_SCORE_MAX = PHRED_64_OFFSET_MAX - PHRED_64_OFFSET_MIN;
52
53	//! Offset used by Solexa encoding quality scores
54	const int SOLEXA_OFFSET_MIN = '@';
55	//! The maximum ASCII value allowed for encoded Solexa scores
56	const int SOLEXA_OFFSET_MAX = 'h';
57	//! Minimum Phred encoded score allowed
58	const int SOLEXA_SCORE_MIN = -5;
59	//! Maximum Phred encoded score allowed
60	const int SOLEXA_SCORE_MAX = SOLEXA_OFFSET_MAX - SOLEXA_OFFSET_MIN;
61
62	///////////////////////////////////////////////////////////////////////////////
63	// Pre-calculation of Solexa <-> Phred conversions
64
65	std::array<char, 256>
66	calc_solexa_to_phred33()	1✔
67	{
68	std::array<char, 256> scores = { 0 };	1✔
69
70	for (int i = SOLEXA_SCORE_MIN; i <= SOLEXA_SCORE_MAX; ++i) {	47✔
71	const double score = round(10.0 * log10(1.0 + pow(10, (i / 10.0))));	46✔
72	const int transformed =	46✔
73	std::max<int>(PHRED_SCORE_MIN, std::min<int>(PHRED_SCORE_MAX, score));	138✔
74	scores.at(i + SOLEXA_OFFSET_MIN) = PHRED_OFFSET_MIN + transformed;	92✔
75	}
76
77	return scores;	1✔
78	}
79
80	const auto g_solexa_to_phred33 = calc_solexa_to_phred33();
81
82	///////////////////////////////////////////////////////////////////////////////
83
84	/**
85	* Uppercase letters in the range a-z, but mangles other characters. Resulting
86	* values outside of A-Z should be compared against/reported, since they may not
87	* reflect the original value.
88	*/
89	constexpr char
90	mangle_to_upper(const char c)	8,285✔
91	{
92	return c & 0xDF;	8,285✔
93	}
94
95	/** Returns the quality score in a form that is human readable */
96	std::string
97	escape_raw_score(char raw)	20✔
98	{
99	return shell_escape(std::string(1, raw));	80✔
100	}
101
102	[[noreturn]] void
103	throw_invalid_phred_33(const char raw)	1✔
104	{
105	const int score = raw - PHRED_33_OFFSET_MIN;	1✔
106	const int alt_score = raw - PHRED_64_OFFSET_MIN;	1✔
107	const int max_score = PHRED_33_SCORE_MAX;	1✔
108	const auto esc = escape_raw_score(raw);	1✔
109	const auto esc_max = escape_raw_score(PHRED_33_OFFSET_MAX);	1✔
110
111	std::ostringstream ss;	1✔
112	ss << "Found Phred+33 encoded quality score of " << score << " (encoded as "	1✔
113	<< esc << "), which is greater than the expected maximum score of "	2✔
114	<< max_score << " (encoded as " << esc_max << "). This suggests that "	1✔
115	<< "input may be Phred+64 encoded.\n\n"
116
117	<< "If the quality scores are actually Phred+64 encoded, which would "
118	<< "mean that the Phred quality score is " << alt_score	2✔
119	<< ", then use the '--quality-format 64' command-line option.\n\n"
120
121	<< "If the quality scores are Phred+33 encoded, but have higher than "
122	"expected scores, then use '--quality-format sam' to permit the full "
123	"range of quality scores.\n\n"
124
125	<< "See the documentation for more information.";	1✔
126
127	throw fastq_error(ss.str());	4✔
128	}	3✔
129
130	[[noreturn]] void
131	throw_invalid_phred_64(const char raw)	3✔
132	{
133	AR_REQUIRE(raw < SOLEXA_OFFSET_MIN,	3✔
134	"invalid_phred called on valid PHRED score");
135
136	const int score = raw - PHRED_64_OFFSET_MIN;	3✔
137	const int alt_score = raw - PHRED_33_OFFSET_MIN;	3✔
138	const int min_score = PHRED_64_SCORE_MIN;	3✔
139	const auto esc = escape_raw_score(raw);	3✔
140	const auto esc_min = escape_raw_score(PHRED_64_OFFSET_MIN);	3✔
141
142	std::ostringstream ss;	3✔
143	ss << "Found Phred+64 encoded quality score of " << score << " (encoded as "	3✔
144	<< esc << "), which is less than the expected minimum score of "	6✔
145	<< min_score << " (encoded as " << esc_min << "). This suggests that "	3✔
146	<< "input may be Phred+33 encoded.\n\n"
147
148	<< "If the quality scores are actually Phred+33 encoded, which would "
149	<< "mean that the Phred quality score is " << alt_score	6✔
150	<< ", then omit the '--quality-format' command-line option or use "
151	<< "'--quality-format 33' to explicitly set the format.\n\n";	3✔
152
153	if (score >= SOLEXA_SCORE_MIN) {	3✔
154	ss << "The quality score could also be the older Solexa format, which has "	2✔
155	<< "a minimum score of -5, but data of this type is rare. If it is "
156	<< "actually Solexa encoded FASTQ data, then use the '--quality-format "
157	<< "solexa' command-line option.\n\n";	2✔
158	}
159
160	ss << "See the documentation for more information.";	3✔
161
162	throw fastq_error(ss.str());	12✔
163	}	9✔
164
165	[[noreturn]] void
166	throw_invalid_solexa(const char raw)	2✔
167	{
168	const int score = raw - SOLEXA_OFFSET_MIN;	2✔
169	const int alt_score = raw - PHRED_33_OFFSET_MIN;	2✔
170	const int min_score = SOLEXA_SCORE_MIN;	2✔
171	const auto esc = escape_raw_score(raw);	2✔
172	const auto esc_min = escape_raw_score(SOLEXA_OFFSET_MIN + SOLEXA_SCORE_MIN);	2✔
173	const auto esc_max = escape_raw_score(SOLEXA_OFFSET_MAX);	2✔
174
175	std::ostringstream ss;	2✔
176	if (score < SOLEXA_SCORE_MIN) {	2✔
177	ss << "Found Solexa encoded quality score of " << score << " (encoded as "	1✔
178	<< esc << "), which is less than the expected minimum score of "	2✔
179	<< min_score << " (encoded as " << esc_min << "). This suggests that "	1✔
180	<< "input may be Phred+33 encoded.\n\n"
181
182	<< "If the quality scores are actually Phred+33 encoded, which would "
183	<< "mean that the Phred quality score is " << alt_score	2✔
184	<< ", then omit the '--quality-format' command-line option or use "
185	<< "'--quality-format 33' to explicitly set the format.\n\n"
186
187	<< "See the documentation for more information.";	1✔
188	} else if (raw > SOLEXA_OFFSET_MAX) {	1✔
189	ss << "Found Solexa encoded quality score of " << score << " (encoded as "	1✔
190	<< esc << "), which is greater than the expected maximum score of "	2✔
191	<< SOLEXA_SCORE_MAX << " (encoded as " << esc_max << ").\n\n"	1✔
192
193	<< "See the documentation for more information.";	2✔
194	} else {
195	AR_FAIL("invalid_phred called on valid PHRED score");	×
196	}
197
198	throw fastq_error(ss.str());	8✔
199	}	8✔
200
201	[[noreturn]] void
202	throw_invalid_score(const quality_encoding encoding, const char raw_score)	12✔
203	{
204	if (raw_score < PHRED_OFFSET_MIN \|\| raw_score > PHRED_OFFSET_MAX) {	12✔
205	std::ostringstream ss;	6✔
206
207	ss << "Found raw FASTQ quality score of " << static_cast<int>(raw_score)	6✔
208	<< " (" << escape_raw_score(raw_score) << "). This is outside of the "	12✔
209	<< "range of valid ASCII encoded quality scores (" << PHRED_OFFSET_MIN	12✔
210	<< " to " << PHRED_OFFSET_MAX << "), meaning that the input file is "	6✔
211	<< "either corrupt or not in FASTQ format!";	6✔
212
213	throw fastq_error(ss.str());	24✔
214	}	6✔
215
216	switch (encoding) {	6✔
217	case quality_encoding::phred_33:	1✔
218	throw_invalid_phred_33(raw_score);	1✔
219
220	case quality_encoding::phred_64:	3✔
221	throw_invalid_phred_64(raw_score);	3✔
222
223	case quality_encoding::solexa:	2✔
224	throw_invalid_solexa(raw_score);	2✔
225
226	case quality_encoding::sam:	×
227	AR_FAIL("This case should have been handled by initial check");	×
228
229	default:	×
230	AR_FAIL("Invalid quality encoding");	×
231	}
232	}
233
234	[[noreturn]] void
235	throw_invalid_base(char c)	4✔
236	{
237	std::ostringstream stream;	4✔
238
239	switch (mangle_to_upper(c)) {	4✔
240	case 'B': // C / G / T	2✔
241	case 'D': // A / G / T	2✔
242	case 'H': // A / C / T	2✔
243	case 'K': // G / T	2✔
244	case 'M': // A / C	2✔
245	case 'R': // A / G	2✔
246	case 'S': // C / G	2✔
247	case 'V': // A / C / G	2✔
248	case 'W': // A / T	2✔
249	case 'Y': // C / T	2✔
250	stream << "found degenerate base '" << c << "' in FASTQ sequence, but "	2✔
251	<< "only bases A, C, G, T and N are supported. Use the option "
252	"--mask-degenerate-bases to convert degenerate bases to N";	2✔
253	break;
254
255	case 'U': // Uracils	×
256	stream << "found uracil (U) in FASTQ sequence, but only bases A, C, G, "	×
257	"T and N are supported. Use the option --convert-uracils to "
258	"convert uracils (U) to thymine (T)";	×
259	break;
260
261	default:	2✔
262	stream << "invalid character " << log_escape(std::string(1, c))	10✔
263	<< " found in FASTQ sequence";	4✔
264	break;	2✔
265	}
266
267	throw fastq_error(stream.str());	16✔
268	}	4✔
269
270	void
271	process_nucleotides_strict(std::string& nucleotides)	1,333✔
272	{
273	for (char& nuc : nucleotides) {	30,167✔
274	// Fast ASCII letter uppercase
275	const auto upper_case = mangle_to_upper(nuc);	8,281✔
276	switch (upper_case) {	8,281✔
277	case 'A':	8,277✔
278	case 'C':	8,277✔
279	case 'G':	8,277✔
280	case 'T':	8,277✔
281	case 'N':	8,277✔
282	nuc = upper_case;	8,277✔
283	break;	8,277✔
284
285	default:	4✔
286	throw_invalid_base(nuc);	4✔
287	}
288	}
289	}	1,329✔
290
291	void
292	process_nucleotides_lenient(std::string& nucleotides,	×
293	const bool convert_uracil,
294	const bool mask_degenerate)
295	{
296	for (char& nuc : nucleotides) {	×
297	// Fast ASCII letter uppercase
NEW 298	const auto upper_case = mangle_to_upper(nuc);	×
299	switch (upper_case) {	×
300	case 'A':	×
301	case 'C':	×
302	case 'G':	×
303	case 'T':	×
304	case 'N':	×
305	nuc = upper_case;	×
306	break;	×
307
308	// Uracils
309	case 'U':	×
310	if (convert_uracil) {	×
311	nuc = 'T';	×
312	break;	×
313	} else {
314	throw_invalid_base(nuc);	×
315	}
316
317	// IUPAC encoded degenerate bases
318	case 'B': // C / G / T	×
319	case 'D': // A / G / T	×
320	case 'H': // A / C / T	×
321	case 'K': // G / T	×
322	case 'M': // A / C	×
323	case 'R': // A / G	×
324	case 'S': // C / G	×
325	case 'V': // A / C / G	×
326	case 'W': // A / T	×
327	case 'Y': // C / T	×
328	if (mask_degenerate) {	×
329	nuc = 'N';	×
330	break;	×
331	} else {
332	throw_invalid_base(nuc);	×
333	}
334
335	default:	×
336	throw_invalid_base(nuc);	×
337	}
338	}
339	}
340
341	} // namespace
342
343	///////////////////////////////////////////////////////////////////////////////
344
345	fastq_encoding::fastq_encoding(quality_encoding encoding,	48✔
346	degenerate_encoding degenerate,
347	uracil_encoding uracils) noexcept	48✔
348	: m_mask_degenerate()	48✔
349	, m_convert_uracil()	48✔
350	, m_encoding(encoding)	48✔
351	, m_offset_min()	48✔
352	, m_offset_max()	48✔
353	{
354	switch (degenerate) {	48✔
355	case degenerate_encoding::reject:	48✔
356	m_mask_degenerate = false;	48✔
357	break;	48✔
358	case degenerate_encoding::mask:	×
359	m_mask_degenerate = true;	×
360	break;	×
361	default:	×
362	AR_FAIL("invalid degenerate encoding value");	×
363	}
364
365	switch (uracils) {	48✔
366	case uracil_encoding::convert:	×
367	m_convert_uracil = true;	×
368	break;	×
369	case uracil_encoding::reject:	48✔
370	m_convert_uracil = false;	48✔
371	break;	48✔
372	default:	×
373	AR_FAIL("invalid uracil encoding value");	×
374	}
375
376	switch (encoding) {	48✔
377	case quality_encoding::phred_33:	12✔
378	m_offset_min = PHRED_33_OFFSET_MIN;	12✔
379	m_offset_max = PHRED_33_OFFSET_MAX;	12✔
380	break;	12✔
381
382	case quality_encoding::phred_64:	12✔
383	m_offset_min = PHRED_64_OFFSET_MIN;	12✔
384	m_offset_max = PHRED_64_OFFSET_MAX;	12✔
385	break;	12✔
386
387	case quality_encoding::solexa:	12✔
388	m_offset_min = SOLEXA_OFFSET_MIN;	12✔
389	m_offset_max = SOLEXA_OFFSET_MAX;	12✔
390	break;	12✔
391
392	case quality_encoding::sam:	12✔
393	m_offset_min = PHRED_OFFSET_MIN;	12✔
394	m_offset_max = PHRED_OFFSET_MAX;	12✔
395	break;	12✔
396
397	default:	×
398	AR_FAIL("unknown quality encoding");	×
399	}
400	}	48✔
401
402	void
403	fastq_encoding::process_nucleotides(std::string& sequence) const	1,333✔
404	{
405	if (m_mask_degenerate \|\| m_convert_uracil) {	1,333✔
406	process_nucleotides_lenient(sequence, m_convert_uracil, m_mask_degenerate);	×
407	} else {
408	process_nucleotides_strict(sequence);	1,333✔
409	}
410	}	1,329✔
411
412	void
413	fastq_encoding::process_qualities(std::string& qualities) const	949✔
414	{
415	if (m_encoding == quality_encoding::solexa) {	949✔
416	for (auto& quality : qualities) {	94✔
417	const char current = quality;	20✔
418	// TODO: Handle negative values, e.g. è
419	if (!(quality = g_solexa_to_phred33.at(quality))) {	40✔
420	throw_invalid_score(m_encoding, current);	2✔
421	}
422	}
423	} else {
424	for (auto& quality : qualities) {	23,277✔
425	if (quality < m_offset_min \|\| quality > m_offset_max) {	6,507✔
426	throw_invalid_score(m_encoding, quality);	10✔
427	}
428
429	// Convert Phred+64 to Phred+33 if needed
430	quality -= m_offset_min - PHRED_33_OFFSET_MIN;	6,497✔
431	}
432	}
433	}	937✔
434
435	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #72

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