#110

Committed 30 Apr 2025 09:26AM UTC coverage: 66.988% (-0.03%) from 67.02%

Build # #110

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

improve estimation of gzip buffer capacity needs (#136)

This reduces the amount of allocations and fixes compression failures
when using --compression-level 1

Coverage Stats

0 of 14 new or added lines in 1 file covered. (0.0%)

9742 of 14543 relevant lines covered (66.99%)

3042.04 hits per line

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

0.0

/src/fastq_io.cpp

// SPDX-License-Identifier: GPL-3.0-or-later
// SPDX-FileCopyrightText: 2015 Mikkel Schubert <mikkelsch@gmail.com>
#include "fastq_io.hpp"      // declarations
#include "debug.hpp"         // for AR_REQUIRE, AR_REQUIRE_SINGLE_THREAD
#include "errors.hpp"        // for io_error, gzip_error, fastq_error
#include "fastq.hpp"         // for fastq
#include "fastq_enc.hpp"     // for MATE_SEPARATOR
#include "output.hpp"        // for output_file
#include "simd.hpp"          // for size_t
#include "statistics.hpp"    // for fastq_statistics, fastq_stats_ptr, stat...
#include "strutils.hpp"      // for shell_escape, string_vec, ends_with
#include "userconfig.hpp"    // for userconfig
#include "utilities.hpp"     // for prng_seed
#include <algorithm>         // for max, min
#include <cerrno>            // for errno
#include <cstring>           // for size_t, memcpy
#include <isa-l/crc.h>       // for crc32_gzip_refl
#include <isa-l/igzip_lib.h> // for isal_zstream, isal_deflate_init, isal_d...
#include <libdeflate.h>      // for libdeflate_alloc_compressor, libdeflate...
#include <memory>            // for unique_ptr, make_unique, __shared_ptr_a...
#include <sstream>           // for basic_ostream, basic_ostringstream, ope...
#include <utility>           // for move, swap

namespace adapterremoval {

// Default bgzip header, as described in the SAM spec. v1.6 section 4.1.
// Includes 2 trailing placeholder bytes for total block size (BSIZE)
constexpr std::string_view BGZF_HEADER = {
  "\37\213\10\4\0\0\0\0\0\377\6\0\102\103\2\0\0\0",
  18
};

// Eof of file marker for bgzip files; see SAM spec. v1.6 section 4.1.2
constexpr std::string_view BGZF_EOF = {
  "\37\213\10\4\0\0\0\0\0\377\6\0\102\103\2\0\33\0\3\0\0\0\0\0\0\0\0\0",
  28,
};

////////////////////////////////////////////////////////////////////////////////
// Helper function for isa-l

namespace {

//! The compression level used for block/stream compression with isa-l
constexpr size_t ISAL_COMPRESSION_LEVEL = 1;
//! The default buffer size for compression at level ISAL_COMPRESSION_LEVEL
constexpr size_t ISAL_BUFFER_SIZE = ISAL_DEF_LVL1_SMALL;

/**
 * ISA-l streaming is enabled only at compression level 1, since little
 * difference was observed between levels 1 to 3. However, it still offers a
 * faster compression (with a lower ratio) than libdeflate level 1.
 */
bool
is_isal_streaming_enabled(const output_file file, unsigned compression_level)
{
  switch (file.format) {
    case output_format::fastq:
    case output_format::sam:
    case output_format::bam:
    case output_format::ubam:
      return false;
    case output_format::fastq_gzip:
    case output_format::sam_gzip:
      return compression_level == ISAL_COMPRESSION_LEVEL;
    default:
      AR_FAIL("invalid output format");
  }
}

} // namespace

///////////////////////////////////////////////////////////////////////////////
// Implementations for 'read_fastq'

enum class read_fastq::file_type
{
  read_1,
  read_2,
  interleaved
};

namespace {

bool
read_record(joined_line_readers& reader, fastq_vec& chunk)
{
  // Line numbers change as we attempt to read the record, and potentially
  // points to the next record in the case of invalid qualities/nucleotides
  const auto line_number = reader.linenumber();

  try {
    chunk.emplace_back();
    auto& record = chunk.back();

    if (record.read_unsafe(reader)) {
      return true;
    } else {
      chunk.pop_back();
      return false;
    }
  } catch (const fastq_error& error) {
    std::ostringstream stream;
    stream << "Error reading FASTQ record from '" << reader.filename()
           << "' at line " << line_number << "; aborting:\n"
           << indent_lines(error.what());

    throw fastq_error(stream.str());
  }
}

const string_vec&
select_filenames(const userconfig& config, const read_fastq::file_type mode)
{
  switch (mode) {
    case read_fastq::file_type::read_1:
    case read_fastq::file_type::interleaved:
      return config.input_files_1;
    case read_fastq::file_type::read_2:
      return config.input_files_2;
    default:
      AR_FAIL("invalid read_fastq::file_type value");
  }
}

/** Estimates an upper bound for the required capacity for gzip compression */
constexpr size_t
estimate_capacity(size_t input_size, bool eof)
{
  return input_size                    //
         + BGZF_HEADER.size()          // Standard bgzip header
         + 1                           // BFINAL | BTYPE
         + 4                           // LEN + NLEN
         + 4                           // CRC32
         + 4                           // ISIZE
         + (eof ? BGZF_EOF.size() : 0) // Standard bgzip tail
    ;
}

} // namespace

read_fastq::read_fastq(const userconfig& config,
                       const size_t next_step,
                       const read_fastq::file_type mode,
                       statistics& stats)
  : analytical_step(processing_order::ordered, "read_fastq")
  , m_reader(select_filenames(config, mode))
  , m_next_step(next_step)
  , m_mode(mode)
  , m_head(config.head)
  , m_mate_separator(config.mate_separator)
  , m_mate_separator_identified(config.mate_separator)
  , m_duplication_1(stats.duplication_1)
  , m_duplication_2(stats.duplication_2)
{
  AR_REQUIRE(m_duplication_1 && m_duplication_2);
}

void
read_fastq::add_steps(scheduler& sch,
                      const userconfig& config,
                      size_t next_step,
                      statistics& stats)
{
  const auto add_step = [&sch, &config, &stats](auto next_step, auto type) {
    return sch.add<read_fastq>(config, next_step, type, stats);
  };

  if (config.interleaved_input) {
    add_step(next_step, read_fastq::file_type::interleaved);
  } else if (config.paired_ended_mode) {
    next_step = add_step(next_step, read_fastq::file_type::read_2);
    add_step(next_step, read_fastq::file_type::read_1);
  } else {
    add_step(next_step, read_fastq::file_type::read_1);
  }
}

chunk_vec
read_fastq::process(chunk_ptr chunk)
{
  AR_REQUIRE_SINGLE_THREAD(m_lock);

  if (m_mode == file_type::read_1 || m_mode == file_type::interleaved) {
    // The scheduler only terminates when the first step stops returning chunks
    if (m_eof) {
      return {};
    }

    chunk = std::make_unique<analytical_chunk>();
  } else {
    AR_REQUIRE(!m_eof && chunk);
  }

  auto& reads_1 = chunk->reads_1;
  auto& reads_2 = chunk->reads_2;

  if (m_mode == file_type::read_1 || m_mode == file_type::interleaved) {
    if (m_mode == file_type::read_1) {
      read_single_end(reads_1, *m_duplication_1);
    } else {
      read_interleaved(reads_1, reads_2);
    }
  } else if (m_mode == file_type::read_2) {
    read_single_end(reads_2, *m_duplication_2);
  } else {
    AR_FAIL("invalid file_type value");
  }

  if (m_mode != file_type::read_1) {
    if (reads_1.size() != reads_2.size()) {
      throw fastq_error("Found unequal number of mate 1 and mate 2 reads; "
                        "input files may be truncated. Please fix before "
                        "continuing.");
    } else if (!m_mate_separator_identified) {
      AR_REQUIRE(reads_1.size() == reads_2.size());
      // Mate separators are identified using the first block, in order to
      // reduce the need for locking in the post-processing step

      // Attempt to determine the mate separator character
      m_mate_separator = fastq::guess_mate_separator(reads_1, reads_2);
      m_mate_separator_identified = true;
    }
  }

  // Head must be checked after the first loop, to produce at least one chunk
  m_eof |= !m_head;
  chunk->eof = m_eof;
  chunk->mate_separator = m_mate_separator;
  chunk->first = m_first;
  m_first = false;

  chunk_vec chunks;
  chunks.emplace_back(m_next_step, std::move(chunk));
  return chunks;
}

void
read_fastq::read_single_end(fastq_vec& reads, duplication_statistics& stats)
{
  for (; reads.size() < INPUT_READS && m_head && !m_eof; m_head--) {
    if (read_record(m_reader, reads)) {
      stats.process(reads.back());
    } else {
      m_eof = true;
    }
  }
}

void
read_fastq::read_interleaved(fastq_vec& reads_1, fastq_vec& reads_2)
{
  for (; reads_1.size() < INPUT_READS && m_head && !m_eof; m_head--) {
    if (read_record(m_reader, reads_1)) {
      m_duplication_1->process(reads_1.back());
    } else {
      m_eof = true;
      break;
    }

    if (read_record(m_reader, reads_2)) {
      m_duplication_2->process(reads_2.back());
    }
  }
}

void
read_fastq::finalize()
{
  AR_REQUIRE_SINGLE_THREAD(m_lock);
  AR_REQUIRE(m_eof);
}

///////////////////////////////////////////////////////////////////////////////
// Implementations for 'post_process_fastq'

post_process_fastq::post_process_fastq(const userconfig& config,
                                       size_t next_step,
                                       statistics& stats)
  : analytical_step(processing_order::unordered, "post_process_fastq")
  , m_statistics_1(stats.input_1)
  , m_statistics_2(stats.input_2)
  , m_next_step(next_step)
  , m_encoding(config.io_encoding)
  , m_timer(config.log_progress)
{
  AR_REQUIRE(m_statistics_1 && m_statistics_2);

  for (size_t i = 0; i < config.max_threads; ++i) {
    m_stats.emplace_back(config.report_sample_rate, prng_seed());
  }
}

chunk_vec
post_process_fastq::process(chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  auto& reads_1 = chunk->reads_1;
  auto& reads_2 = chunk->reads_2;

  auto stats = m_stats.acquire();

  AR_REQUIRE((reads_1.size() == reads_2.size()) || reads_2.empty());
  if (reads_2.empty()) {
    for (auto& read_1 : reads_1) {
      read_1.post_process(m_encoding);
      stats->stats_1.process(read_1);
    }
  } else {
    auto it_1 = reads_1.begin();
    auto it_2 = reads_2.begin();
    for (; it_1 != reads_1.end(); ++it_1, ++it_2) {
      fastq::normalize_paired_reads(*it_1, *it_2, chunk->mate_separator);

      it_1->post_process(m_encoding);
      stats->stats_1.process(*it_1);

      it_2->post_process(m_encoding);
      stats->stats_2.process(*it_2);
    }

    // fastq::normalize_paired_reads replaces the mate separator if present
    if (chunk->mate_separator) {
      chunk->mate_separator = MATE_SEPARATOR;
    }
  }

  m_stats.release(stats);

  {
    std::unique_lock<std::mutex> lock(m_timer_lock);
    m_timer.increment(reads_1.size() + reads_2.size());
  }

  chunk_vec chunks;
  chunks.emplace_back(m_next_step, std::move(chunk));

  return chunks;
}

void
post_process_fastq::finalize()
{
  AR_REQUIRE_SINGLE_THREAD(m_timer_lock);

  while (auto it = m_stats.try_acquire()) {
    *m_statistics_1 += it->stats_1;
    *m_statistics_2 += it->stats_2;
  }

  m_timer.finalize();
}

///////////////////////////////////////////////////////////////////////////////
// Implementations for 'split_fastq'

split_fastq::split_fastq(const userconfig& config,
                         const output_file& file,
                         size_t next_step)
  : analytical_step(processing_order::ordered, "split_fastq")
  , m_next_step(next_step)
  , m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))
{
}

void
split_fastq::finalize()
{
  AR_REQUIRE_SINGLE_THREAD(m_lock);

  AR_REQUIRE(m_eof);
  AR_REQUIRE(m_buffer.capacity() == 0);
}

chunk_vec
split_fastq::process(const chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  AR_REQUIRE_SINGLE_THREAD(m_lock);
  AR_REQUIRE(!m_eof);
  m_eof = chunk->eof;

  chunk_vec chunks;
  for (const auto& src : chunk->buffers) {
    for (size_t src_offset = 0; src_offset < src.size();) {
      const auto n =
        std::min(src.size() - src_offset, BGZF_BLOCK_SIZE - m_buffer.size());
      m_buffer.append(src.data() + src_offset, n);

      src_offset += n;

      if (m_buffer.size() == BGZF_BLOCK_SIZE) {
        auto block = std::make_unique<analytical_chunk>();

        if (m_isal_stream) {
          m_isal_crc32 =
            crc32_gzip_refl(m_isal_crc32, m_buffer.data(), m_buffer.size());
        }

        block->uncompressed_size = m_buffer.size();
        block->buffers.emplace_back(std::move(m_buffer));

        chunks.emplace_back(m_next_step, std::move(block));

        m_buffer = buffer();
        m_buffer.reserve(BGZF_BLOCK_SIZE);
      }
    }
  }

  if (m_eof) {
    auto block = std::make_unique<analytical_chunk>();
    block->eof = true;

    if (m_isal_stream) {
      m_isal_crc32 =
        crc32_gzip_refl(m_isal_crc32, m_buffer.data(), m_buffer.size());
    }

    block->crc32 = m_isal_crc32;
    block->uncompressed_size = m_buffer.size();
    block->buffers.emplace_back(std::move(m_buffer));

    chunks.emplace_back(m_next_step, std::move(block));
  }

  return chunks;
}

///////////////////////////////////////////////////////////////////////////////
// Implementations for 'gzip_split_fastq'

namespace {

size_t
isal_deflate_block(buffer& input_buffer,
                   buffer& output_buffer,
                   const size_t output_offset,
                   const bool eof)
{
  isal_zstream stream{};
  isal_deflate_stateless_init(&stream);

  stream.flush = FULL_FLUSH;
  stream.end_of_stream = eof;

  stream.level = ISAL_COMPRESSION_LEVEL;
  stream.level_buf_size = ISAL_BUFFER_SIZE;
  buffer level_buffer{ stream.level_buf_size };
  stream.level_buf = level_buffer.data();

  stream.avail_in = input_buffer.size();
  stream.next_in = input_buffer.data();
  stream.next_out = output_buffer.data() + output_offset;
  stream.avail_out = output_buffer.size() - output_offset;

  const auto ec = isal_deflate_stateless(&stream);
  switch (ec) {
    case COMP_OK:
      break;
    case INVALID_FLUSH:
      throw gzip_error("isal_deflate_stateless: invalid flush");
    case ISAL_INVALID_LEVEL:
      throw gzip_error("isal_deflate_stateless: invalid level");
    case ISAL_INVALID_LEVEL_BUF:
      throw gzip_error("isal_deflate_stateless: invalid buffer size");
    default: {
      std::ostringstream os;
      os << "isal_deflate_stateless: unknown error " << ec;

      throw gzip_error(os.str());
    }
  }

  // The easily compressible input should fit in a single output block
  AR_REQUIRE(stream.avail_in == 0);

  return stream.total_out;
}

} // namespace

gzip_split_fastq::gzip_split_fastq(const userconfig& config,
                                   const output_file& file,
                                   size_t next_step)
  : analytical_step(processing_order::unordered, "gzip_split_fastq")
  , m_config(config)
  , m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))
  , m_format(file.format)
  , m_next_step(next_step)
{
}

chunk_vec
gzip_split_fastq::process(chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  AR_REQUIRE(chunk->buffers.size() == 1);

  buffer& input_buffer = chunk->buffers.front();
  buffer output_buffer;

  if (m_isal_stream) {
    output_buffer.resize(estimate_capacity(input_buffer.size(), false));
    const auto output_size =
      isal_deflate_block(input_buffer, output_buffer, 0, chunk->eof);
    output_buffer.resize(output_size);
  } else {
    if (m_format == output_format::ubam || m_config.compression_level == 0) {
      output_buffer.reserve(estimate_capacity(input_buffer.size(), chunk->eof));
      output_buffer.append(BGZF_HEADER);
      output_buffer.append_u8(1); // BFINAL=1, BTYPE=00; see RFC1951
      output_buffer.append_u16(input_buffer.size());
      output_buffer.append_u16(~input_buffer.size());
      output_buffer.append(input_buffer);
    } else if (m_config.compression_level == ISAL_COMPRESSION_LEVEL) {
      output_buffer.reserve(estimate_capacity(input_buffer.size(), chunk->eof));
      output_buffer.append(BGZF_HEADER);
      output_buffer.resize(output_buffer.capacity());

      const auto output_size = isal_deflate_block(input_buffer,
                                                  output_buffer,
                                                  BGZF_HEADER.size(),
                                                  true);

      // Resize the buffer to the actually used size
      output_buffer.resize(output_size + BGZF_HEADER.size());
    } else {
      // Libdeflate compression levels 1 to 12 are mapped onto 2 to 13
      AR_REQUIRE(m_config.compression_level >= 2 &&
                 m_config.compression_level <= 13);
      auto* compressor =
        libdeflate_alloc_compressor(m_config.compression_level - 1);
      const auto output_bound =
        libdeflate_deflate_compress_bound(compressor, input_buffer.size());

      output_buffer.reserve(estimate_capacity(output_bound, chunk->eof));
      output_buffer.append(BGZF_HEADER);
      output_buffer.resize(output_buffer.capacity());

      const auto output_size =
        libdeflate_deflate_compress(compressor,
                                    input_buffer.data(),
                                    input_buffer.size(),
                                    output_buffer.data() + BGZF_HEADER.size(),
                                    output_buffer.size() - BGZF_HEADER.size());
      libdeflate_free_compressor(compressor);
      // The easily compressible input should fit in a single output block
      AR_REQUIRE(output_size);

      // Resize the buffer to the actually used size
      output_buffer.resize(output_size + BGZF_HEADER.size());
    }

    const auto input_crc32 =
      libdeflate_crc32(0, input_buffer.data(), input_buffer.size());
    output_buffer.append_u32(input_crc32);         // checksum of data
    output_buffer.append_u32(input_buffer.size()); // size of data

    AR_REQUIRE(output_buffer.size() <= BGZF_MAX_BLOCK_SIZE);
    // Write the final block size; -1 to fit 65536 in 16 bit
    output_buffer.put_u16(16, output_buffer.size() - 1);

    if (chunk->eof) {
      output_buffer.append(BGZF_EOF);
    }
  }

  // Enable reuse of the analytical_chunks
  std::swap(input_buffer, output_buffer);

  chunk_vec chunks;
  chunks.emplace_back(m_next_step, std::move(chunk));

  return chunks;
}

///////////////////////////////////////////////////////////////////////////////
// Implementations for 'write_fastq'

write_fastq::write_fastq(const userconfig& config, const output_file& file)
  // Allow disk IO and writing to STDOUT at the same time
  : analytical_step(processing_order::ordered_io, "write_fastq")
  , m_output(file.name)
  , m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))
{
  if (m_isal_stream) {
    buffer level_buf(ISAL_BUFFER_SIZE);
    buffer output_buf(OUTPUT_BLOCK_SIZE);

    struct isal_zstream stream = {};
    struct isal_gzip_header header = {};

    isal_gzip_header_init(&header);
    isal_deflate_init(&stream);
    stream.avail_in = 0;
    stream.flush = NO_FLUSH;
    stream.level = ISAL_COMPRESSION_LEVEL;
    stream.level_buf = level_buf.data();
    stream.level_buf_size = level_buf.size();
    stream.gzip_flag = IGZIP_GZIP_NO_HDR;
    stream.next_out = output_buf.data();
    stream.avail_out = output_buf.size();

    const auto ret = isal_write_gzip_header(&stream, &header);
    AR_REQUIRE(ret == 0, "buffer was not large enough for gzip header");

    output_buf.resize(stream.total_out);

    m_output.write(output_buf);
  }
}

chunk_vec
write_fastq::process(chunk_ptr chunk)
{
  AR_REQUIRE(chunk);
  AR_REQUIRE_SINGLE_THREAD(m_lock);
  AR_REQUIRE(!m_eof);

  try {
    m_eof = chunk->eof;
    m_uncompressed_bytes += chunk->uncompressed_size;

    const auto mode = (m_eof && !m_isal_stream) ? flush::on : flush::off;

    m_output.write(chunk->buffers, mode);

    if (m_eof && m_isal_stream) {
      buffer trailer;
      trailer.append_u32(chunk->crc32);
      trailer.append_u32(m_uncompressed_bytes);

      m_output.write(trailer, flush::on);
    }
  } catch (const std::ios_base::failure&) {
    std::ostringstream msg;
    msg << "Error writing to FASTQ file " << shell_escape(m_output.filename());

    throw io_error(msg.str(), errno);
  }

  return {};
}

void
write_fastq::finalize()
{
  AR_REQUIRE_SINGLE_THREAD(m_lock);
  AR_REQUIRE(m_eof);

  // Close file to trigger any exceptions due to badbit / failbit
  try {
    m_output.close();
  } catch (const std::ios_base::failure&) {
    std::ostringstream msg;
    msg << "Error closing FASTQ file " << shell_escape(m_output.filename());

    throw io_error(msg.str(), errno);
  }
}

} // namespace adapterremoval

1	// SPDX-License-Identifier: GPL-3.0-or-later
2	// SPDX-FileCopyrightText: 2015 Mikkel Schubert <mikkelsch@gmail.com>
3	#include "fastq_io.hpp" // declarations
4	#include "debug.hpp" // for AR_REQUIRE, AR_REQUIRE_SINGLE_THREAD
5	#include "errors.hpp" // for io_error, gzip_error, fastq_error
6	#include "fastq.hpp" // for fastq
7	#include "fastq_enc.hpp" // for MATE_SEPARATOR
8	#include "output.hpp" // for output_file
9	#include "simd.hpp" // for size_t
10	#include "statistics.hpp" // for fastq_statistics, fastq_stats_ptr, stat...
11	#include "strutils.hpp" // for shell_escape, string_vec, ends_with
12	#include "userconfig.hpp" // for userconfig
13	#include "utilities.hpp" // for prng_seed
14	#include <algorithm> // for max, min
15	#include <cerrno> // for errno
16	#include <cstring> // for size_t, memcpy
17	#include <isa-l/crc.h> // for crc32_gzip_refl
18	#include <isa-l/igzip_lib.h> // for isal_zstream, isal_deflate_init, isal_d...
19	#include <libdeflate.h> // for libdeflate_alloc_compressor, libdeflate...
20	#include <memory> // for unique_ptr, make_unique, __shared_ptr_a...
21	#include <sstream> // for basic_ostream, basic_ostringstream, ope...
22	#include <utility> // for move, swap
23
24	namespace adapterremoval {
25
26	// Default bgzip header, as described in the SAM spec. v1.6 section 4.1.
27	// Includes 2 trailing placeholder bytes for total block size (BSIZE)
28	constexpr std::string_view BGZF_HEADER = {
29	"\37\213\10\4\0\0\0\0\0\377\6\0\102\103\2\0\0\0",
30	18
31	};
32
33	// Eof of file marker for bgzip files; see SAM spec. v1.6 section 4.1.2
34	constexpr std::string_view BGZF_EOF = {
35	"\37\213\10\4\0\0\0\0\0\377\6\0\102\103\2\0\33\0\3\0\0\0\0\0\0\0\0\0",
36	28,
37	};
38
39	////////////////////////////////////////////////////////////////////////////////
40	// Helper function for isa-l
41
42	namespace {
43
44	//! The compression level used for block/stream compression with isa-l
45	constexpr size_t ISAL_COMPRESSION_LEVEL = 1;
46	//! The default buffer size for compression at level ISAL_COMPRESSION_LEVEL
47	constexpr size_t ISAL_BUFFER_SIZE = ISAL_DEF_LVL1_SMALL;
48
49	/**
50	* ISA-l streaming is enabled only at compression level 1, since little
51	* difference was observed between levels 1 to 3. However, it still offers a
52	* faster compression (with a lower ratio) than libdeflate level 1.
53	*/
54	bool
55	is_isal_streaming_enabled(const output_file file, unsigned compression_level)	×
56	{
57	switch (file.format) {	×
58	case output_format::fastq:
59	case output_format::sam:
60	case output_format::bam:
61	case output_format::ubam:
62	return false;
63	case output_format::fastq_gzip:	×
64	case output_format::sam_gzip:	×
65	return compression_level == ISAL_COMPRESSION_LEVEL;	×
66	default:	×
67	AR_FAIL("invalid output format");	×
68	}
69	}
70
71	} // namespace
72
73	///////////////////////////////////////////////////////////////////////////////
74	// Implementations for 'read_fastq'
75
76	enum class read_fastq::file_type
77	{
78	read_1,
79	read_2,
80	interleaved
81	};
82
83	namespace {
84
85	bool
86	read_record(joined_line_readers& reader, fastq_vec& chunk)	×
87	{
88	// Line numbers change as we attempt to read the record, and potentially
89	// points to the next record in the case of invalid qualities/nucleotides
90	const auto line_number = reader.linenumber();	×
91
92	try {	×
93	chunk.emplace_back();	×
94	auto& record = chunk.back();	×
95
96	if (record.read_unsafe(reader)) {	×
97	return true;
98	} else {
99	chunk.pop_back();	×
100	return false;	×
101	}
102	} catch (const fastq_error& error) {	×
103	std::ostringstream stream;	×
104	stream << "Error reading FASTQ record from '" << reader.filename()	×
105	<< "' at line " << line_number << "; aborting:\n"	×
106	<< indent_lines(error.what());	×
107
108	throw fastq_error(stream.str());	×
109	}	×
110	}
111
112	const string_vec&
113	select_filenames(const userconfig& config, const read_fastq::file_type mode)	×
114	{
115	switch (mode) {	×
116	case read_fastq::file_type::read_1:	×
117	case read_fastq::file_type::interleaved:	×
118	return config.input_files_1;	×
119	case read_fastq::file_type::read_2:	×
120	return config.input_files_2;	×
121	default:	×
122	AR_FAIL("invalid read_fastq::file_type value");	×
123	}
124	}
125
126	/** Estimates an upper bound for the required capacity for gzip compression */
127	constexpr size_t
NEW 128	estimate_capacity(size_t input_size, bool eof)	×
129	{
NEW 130	return input_size //	×
131	+ BGZF_HEADER.size() // Standard bgzip header
132	+ 1 // BFINAL \| BTYPE
133	+ 4 // LEN + NLEN
134	+ 4 // CRC32
135	+ 4 // ISIZE
NEW 136	+ (eof ? BGZF_EOF.size() : 0) // Standard bgzip tail	×
137	;
138	}
139
140	} // namespace
141
142	read_fastq::read_fastq(const userconfig& config,	×
143	const size_t next_step,
144	const read_fastq::file_type mode,
145	statistics& stats)	×
146	: analytical_step(processing_order::ordered, "read_fastq")
147	, m_reader(select_filenames(config, mode))	×
148	, m_next_step(next_step)	×
149	, m_mode(mode)	×
150	, m_head(config.head)	×
151	, m_mate_separator(config.mate_separator)	×
152	, m_mate_separator_identified(config.mate_separator)	×
153	, m_duplication_1(stats.duplication_1)	×
154	, m_duplication_2(stats.duplication_2)	×
155	{
156	AR_REQUIRE(m_duplication_1 && m_duplication_2);	×
157	}
158
159	void
160	read_fastq::add_steps(scheduler& sch,	×
161	const userconfig& config,
162	size_t next_step,
163	statistics& stats)
164	{
165	const auto add_step = [&sch, &config, &stats](auto next_step, auto type) {	×
166	return sch.add<read_fastq>(config, next_step, type, stats);	×
167	};
168
169	if (config.interleaved_input) {	×
170	add_step(next_step, read_fastq::file_type::interleaved);	×
171	} else if (config.paired_ended_mode) {	×
172	next_step = add_step(next_step, read_fastq::file_type::read_2);	×
173	add_step(next_step, read_fastq::file_type::read_1);	×
174	} else {
175	add_step(next_step, read_fastq::file_type::read_1);	×
176	}
177	}
178
179	chunk_vec
180	read_fastq::process(chunk_ptr chunk)	×
181	{
182	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
183
184	if (m_mode == file_type::read_1 \|\| m_mode == file_type::interleaved) {	×
185	// The scheduler only terminates when the first step stops returning chunks
186	if (m_eof) {	×
187	return {};	×
188	}
189
190	chunk = std::make_unique<analytical_chunk>();	×
191	} else {
192	AR_REQUIRE(!m_eof && chunk);	×
193	}
194
195	auto& reads_1 = chunk->reads_1;	×
196	auto& reads_2 = chunk->reads_2;	×
197
198	if (m_mode == file_type::read_1 \|\| m_mode == file_type::interleaved) {	×
199	if (m_mode == file_type::read_1) {	×
200	read_single_end(reads_1, *m_duplication_1);	×
201	} else {
202	read_interleaved(reads_1, reads_2);	×
203	}
204	} else if (m_mode == file_type::read_2) {	×
205	read_single_end(reads_2, *m_duplication_2);	×
206	} else {
207	AR_FAIL("invalid file_type value");	×
208	}
209
210	if (m_mode != file_type::read_1) {	×
211	if (reads_1.size() != reads_2.size()) {	×
212	throw fastq_error("Found unequal number of mate 1 and mate 2 reads; "	×
213	"input files may be truncated. Please fix before "
214	"continuing.");	×
215	} else if (!m_mate_separator_identified) {	×
216	AR_REQUIRE(reads_1.size() == reads_2.size());	×
217	// Mate separators are identified using the first block, in order to
218	// reduce the need for locking in the post-processing step
219
220	// Attempt to determine the mate separator character
221	m_mate_separator = fastq::guess_mate_separator(reads_1, reads_2);	×
222	m_mate_separator_identified = true;	×
223	}
224	}
225
226	// Head must be checked after the first loop, to produce at least one chunk
227	m_eof \|= !m_head;	×
228	chunk->eof = m_eof;	×
229	chunk->mate_separator = m_mate_separator;	×
230	chunk->first = m_first;	×
231	m_first = false;	×
232
233	chunk_vec chunks;	×
234	chunks.emplace_back(m_next_step, std::move(chunk));	×
235	return chunks;	×
236	}
237
238	void
239	read_fastq::read_single_end(fastq_vec& reads, duplication_statistics& stats)	×
240	{
241	for (; reads.size() < INPUT_READS && m_head && !m_eof; m_head--) {	×
242	if (read_record(m_reader, reads)) {	×
243	stats.process(reads.back());	×
244	} else {
245	m_eof = true;	×
246	}
247	}
248	}
249
250	void
251	read_fastq::read_interleaved(fastq_vec& reads_1, fastq_vec& reads_2)	×
252	{
253	for (; reads_1.size() < INPUT_READS && m_head && !m_eof; m_head--) {	×
254	if (read_record(m_reader, reads_1)) {	×
255	m_duplication_1->process(reads_1.back());	×
256	} else {
257	m_eof = true;	×
258	break;	×
259	}
260
261	if (read_record(m_reader, reads_2)) {	×
262	m_duplication_2->process(reads_2.back());	×
263	}
264	}
265	}
266
267	void
268	read_fastq::finalize()	×
269	{
270	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
271	AR_REQUIRE(m_eof);	×
272	}
273
274	///////////////////////////////////////////////////////////////////////////////
275	// Implementations for 'post_process_fastq'
276
277	post_process_fastq::post_process_fastq(const userconfig& config,	×
278	size_t next_step,
279	statistics& stats)	×
280	: analytical_step(processing_order::unordered, "post_process_fastq")
281	, m_statistics_1(stats.input_1)	×
282	, m_statistics_2(stats.input_2)	×
283	, m_next_step(next_step)	×
284	, m_encoding(config.io_encoding)	×
285	, m_timer(config.log_progress)	×
286	{
287	AR_REQUIRE(m_statistics_1 && m_statistics_2);	×
288
289	for (size_t i = 0; i < config.max_threads; ++i) {	×
290	m_stats.emplace_back(config.report_sample_rate, prng_seed());	×
291	}
292	}
293
294	chunk_vec
295	post_process_fastq::process(chunk_ptr chunk)	×
296	{
297	AR_REQUIRE(chunk);	×
298	auto& reads_1 = chunk->reads_1;	×
299	auto& reads_2 = chunk->reads_2;	×
300
301	auto stats = m_stats.acquire();	×
302
303	AR_REQUIRE((reads_1.size() == reads_2.size()) \|\| reads_2.empty());	×
304	if (reads_2.empty()) {	×
305	for (auto& read_1 : reads_1) {	×
306	read_1.post_process(m_encoding);	×
307	stats->stats_1.process(read_1);	×
308	}
309	} else {
310	auto it_1 = reads_1.begin();	×
311	auto it_2 = reads_2.begin();	×
312	for (; it_1 != reads_1.end(); ++it_1, ++it_2) {	×
313	fastq::normalize_paired_reads(it_1, it_2, chunk->mate_separator);	×
314
315	it_1->post_process(m_encoding);	×
316	stats->stats_1.process(*it_1);	×
317
318	it_2->post_process(m_encoding);	×
319	stats->stats_2.process(*it_2);	×
320	}
321
322	// fastq::normalize_paired_reads replaces the mate separator if present
323	if (chunk->mate_separator) {	×
324	chunk->mate_separator = MATE_SEPARATOR;	×
325	}
326	}
327
328	m_stats.release(stats);	×
329
330	{	×
331	std::unique_lock<std::mutex> lock(m_timer_lock);	×
332	m_timer.increment(reads_1.size() + reads_2.size());	×
333	}
334
335	chunk_vec chunks;	×
336	chunks.emplace_back(m_next_step, std::move(chunk));	×
337
338	return chunks;	×
339	}
340
341	void
342	post_process_fastq::finalize()	×
343	{
344	AR_REQUIRE_SINGLE_THREAD(m_timer_lock);	×
345
346	while (auto it = m_stats.try_acquire()) {	×
347	*m_statistics_1 += it->stats_1;	×
348	*m_statistics_2 += it->stats_2;	×
349	}
350
351	m_timer.finalize();	×
352	}
353
354	///////////////////////////////////////////////////////////////////////////////
355	// Implementations for 'split_fastq'
356
357	split_fastq::split_fastq(const userconfig& config,	×
358	const output_file& file,
359	size_t next_step)	×
360	: analytical_step(processing_order::ordered, "split_fastq")
361	, m_next_step(next_step)	×
362	, m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))	×
363	{
364	}
365
366	void
367	split_fastq::finalize()	×
368	{
369	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
370
371	AR_REQUIRE(m_eof);	×
372	AR_REQUIRE(m_buffer.capacity() == 0);	×
373	}
374
375	chunk_vec
376	split_fastq::process(const chunk_ptr chunk)	×
377	{
378	AR_REQUIRE(chunk);	×
379	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
380	AR_REQUIRE(!m_eof);	×
381	m_eof = chunk->eof;	×
382
383	chunk_vec chunks;	×
384	for (const auto& src : chunk->buffers) {	×
385	for (size_t src_offset = 0; src_offset < src.size();) {	×
386	const auto n =	×
387	std::min(src.size() - src_offset, BGZF_BLOCK_SIZE - m_buffer.size());	×
388	m_buffer.append(src.data() + src_offset, n);	×
389
390	src_offset += n;	×
391
392	if (m_buffer.size() == BGZF_BLOCK_SIZE) {	×
393	auto block = std::make_unique<analytical_chunk>();	×
394
395	if (m_isal_stream) {	×
396	m_isal_crc32 =	×
397	crc32_gzip_refl(m_isal_crc32, m_buffer.data(), m_buffer.size());	×
398	}
399
400	block->uncompressed_size = m_buffer.size();	×
401	block->buffers.emplace_back(std::move(m_buffer));	×
402
403	chunks.emplace_back(m_next_step, std::move(block));	×
404
405	m_buffer = buffer();	×
406	m_buffer.reserve(BGZF_BLOCK_SIZE);	×
407	}
408	}
409	}
410
411	if (m_eof) {	×
412	auto block = std::make_unique<analytical_chunk>();	×
413	block->eof = true;	×
414
415	if (m_isal_stream) {	×
416	m_isal_crc32 =	×
417	crc32_gzip_refl(m_isal_crc32, m_buffer.data(), m_buffer.size());	×
418	}
419
420	block->crc32 = m_isal_crc32;	×
421	block->uncompressed_size = m_buffer.size();	×
422	block->buffers.emplace_back(std::move(m_buffer));	×
423
424	chunks.emplace_back(m_next_step, std::move(block));	×
425	}
426
427	return chunks;	×
428	}
429
430	///////////////////////////////////////////////////////////////////////////////
431	// Implementations for 'gzip_split_fastq'
432
433	namespace {
434
435	size_t
436	isal_deflate_block(buffer& input_buffer,	×
437	buffer& output_buffer,
438	const size_t output_offset,
439	const bool eof)
440	{
441	isal_zstream stream{};	×
442	isal_deflate_stateless_init(&stream);	×
443
444	stream.flush = FULL_FLUSH;	×
445	stream.end_of_stream = eof;	×
446
447	stream.level = ISAL_COMPRESSION_LEVEL;	×
448	stream.level_buf_size = ISAL_BUFFER_SIZE;	×
449	buffer level_buffer{ stream.level_buf_size };	×
450	stream.level_buf = level_buffer.data();	×
451
452	stream.avail_in = input_buffer.size();	×
453	stream.next_in = input_buffer.data();	×
454	stream.next_out = output_buffer.data() + output_offset;	×
455	stream.avail_out = output_buffer.size() - output_offset;	×
456
NEW 457	const auto ec = isal_deflate_stateless(&stream);	×
NEW 458	switch (ec) {	×
459	case COMP_OK:	×
460	break;	×
461	case INVALID_FLUSH:	×
462	throw gzip_error("isal_deflate_stateless: invalid flush");	×
463	case ISAL_INVALID_LEVEL:	×
464	throw gzip_error("isal_deflate_stateless: invalid level");	×
465	case ISAL_INVALID_LEVEL_BUF:	×
466	throw gzip_error("isal_deflate_stateless: invalid buffer size");	×
NEW 467	default: {	×
NEW 468	std::ostringstream os;	×
NEW 469	os << "isal_deflate_stateless: unknown error " << ec;	×
470
NEW 471	throw gzip_error(os.str());	×
NEW 472	}	×
473	}
474
475	// The easily compressible input should fit in a single output block
476	AR_REQUIRE(stream.avail_in == 0);	×
477
478	return stream.total_out;	×
479	}
480
481	} // namespace
482
483	gzip_split_fastq::gzip_split_fastq(const userconfig& config,	×
484	const output_file& file,
485	size_t next_step)	×
486	: analytical_step(processing_order::unordered, "gzip_split_fastq")
487	, m_config(config)	×
488	, m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))	×
489	, m_format(file.format)	×
490	, m_next_step(next_step)	×
491	{
492	}
493
494	chunk_vec
495	gzip_split_fastq::process(chunk_ptr chunk)	×
496	{
497	AR_REQUIRE(chunk);	×
498	AR_REQUIRE(chunk->buffers.size() == 1);	×
499
500	buffer& input_buffer = chunk->buffers.front();	×
501	buffer output_buffer;	×
502
503	if (m_isal_stream) {	×
NEW 504	output_buffer.resize(estimate_capacity(input_buffer.size(), false));	×
505	const auto output_size =	×
506	isal_deflate_block(input_buffer, output_buffer, 0, chunk->eof);	×
507	output_buffer.resize(output_size);	×
508	} else {
509	if (m_format == output_format::ubam \|\| m_config.compression_level == 0) {	×
NEW 510	output_buffer.reserve(estimate_capacity(input_buffer.size(), chunk->eof));	×
511	output_buffer.append(BGZF_HEADER);	×
512	output_buffer.append_u8(1); // BFINAL=1, BTYPE=00; see RFC1951	×
513	output_buffer.append_u16(input_buffer.size());	×
514	output_buffer.append_u16(~input_buffer.size());	×
515	output_buffer.append(input_buffer);	×
516	} else if (m_config.compression_level == ISAL_COMPRESSION_LEVEL) {	×
NEW 517	output_buffer.reserve(estimate_capacity(input_buffer.size(), chunk->eof));	×
518	output_buffer.append(BGZF_HEADER);	×
519	output_buffer.resize(output_buffer.capacity());	×
520
521	const auto output_size = isal_deflate_block(input_buffer,	×
522	output_buffer,
523	BGZF_HEADER.size(),
524	true);
525
526	// Resize the buffer to the actually used size
527	output_buffer.resize(output_size + BGZF_HEADER.size());	×
528	} else {
529	// Libdeflate compression levels 1 to 12 are mapped onto 2 to 13
530	AR_REQUIRE(m_config.compression_level >= 2 &&	×
531	m_config.compression_level <= 13);
532	auto* compressor =	×
533	libdeflate_alloc_compressor(m_config.compression_level - 1);	×
534	const auto output_bound =	×
535	libdeflate_deflate_compress_bound(compressor, input_buffer.size());	×
536
NEW 537	output_buffer.reserve(estimate_capacity(output_bound, chunk->eof));	×
538	output_buffer.append(BGZF_HEADER);	×
539	output_buffer.resize(output_buffer.capacity());	×
540
541	const auto output_size =	×
542	libdeflate_deflate_compress(compressor,	×
543	input_buffer.data(),	×
544	input_buffer.size(),	×
545	output_buffer.data() + BGZF_HEADER.size(),	×
546	output_buffer.size() - BGZF_HEADER.size());	×
547	libdeflate_free_compressor(compressor);	×
548	// The easily compressible input should fit in a single output block
549	AR_REQUIRE(output_size);	×
550
551	// Resize the buffer to the actually used size
552	output_buffer.resize(output_size + BGZF_HEADER.size());	×
553	}
554
555	const auto input_crc32 =	×
556	libdeflate_crc32(0, input_buffer.data(), input_buffer.size());	×
557	output_buffer.append_u32(input_crc32); // checksum of data	×
558	output_buffer.append_u32(input_buffer.size()); // size of data	×
559
560	AR_REQUIRE(output_buffer.size() <= BGZF_MAX_BLOCK_SIZE);	×
561	// Write the final block size; -1 to fit 65536 in 16 bit
562	output_buffer.put_u16(16, output_buffer.size() - 1);	×
563
564	if (chunk->eof) {	×
565	output_buffer.append(BGZF_EOF);	×
566	}
567	}
568
569	// Enable reuse of the analytical_chunks
570	std::swap(input_buffer, output_buffer);	×
571
572	chunk_vec chunks;	×
573	chunks.emplace_back(m_next_step, std::move(chunk));	×
574
575	return chunks;	×
576	}
577
578	///////////////////////////////////////////////////////////////////////////////
579	// Implementations for 'write_fastq'
580
581	write_fastq::write_fastq(const userconfig& config, const output_file& file)	×
582	// Allow disk IO and writing to STDOUT at the same time
583	: analytical_step(processing_order::ordered_io, "write_fastq")
584	, m_output(file.name)	×
585	, m_isal_stream(is_isal_streaming_enabled(file, config.compression_level))	×
586	{
587	if (m_isal_stream) {	×
588	buffer level_buf(ISAL_BUFFER_SIZE);	×
589	buffer output_buf(OUTPUT_BLOCK_SIZE);	×
590
591	struct isal_zstream stream = {};	×
592	struct isal_gzip_header header = {};	×
593
594	isal_gzip_header_init(&header);	×
595	isal_deflate_init(&stream);	×
596	stream.avail_in = 0;	×
597	stream.flush = NO_FLUSH;	×
598	stream.level = ISAL_COMPRESSION_LEVEL;	×
599	stream.level_buf = level_buf.data();	×
600	stream.level_buf_size = level_buf.size();	×
601	stream.gzip_flag = IGZIP_GZIP_NO_HDR;	×
602	stream.next_out = output_buf.data();	×
603	stream.avail_out = output_buf.size();	×
604
605	const auto ret = isal_write_gzip_header(&stream, &header);	×
606	AR_REQUIRE(ret == 0, "buffer was not large enough for gzip header");	×
607
608	output_buf.resize(stream.total_out);	×
609
610	m_output.write(output_buf);	×
611	}
612	}
613
614	chunk_vec
615	write_fastq::process(chunk_ptr chunk)	×
616	{
617	AR_REQUIRE(chunk);	×
618	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
619	AR_REQUIRE(!m_eof);	×
620
621	try {	×
622	m_eof = chunk->eof;	×
623	m_uncompressed_bytes += chunk->uncompressed_size;	×
624
625	const auto mode = (m_eof && !m_isal_stream) ? flush::on : flush::off;	×
626
627	m_output.write(chunk->buffers, mode);	×
628
629	if (m_eof && m_isal_stream) {	×
630	buffer trailer;	×
631	trailer.append_u32(chunk->crc32);	×
632	trailer.append_u32(m_uncompressed_bytes);	×
633
634	m_output.write(trailer, flush::on);	×
635	}
636	} catch (const std::ios_base::failure&) {	×
637	std::ostringstream msg;	×
638	msg << "Error writing to FASTQ file " << shell_escape(m_output.filename());	×
639
640	throw io_error(msg.str(), errno);	×
641	}	×
642
643	return {};	×
644	}
645
646	void
647	write_fastq::finalize()	×
648	{
649	AR_REQUIRE_SINGLE_THREAD(m_lock);	×
650	AR_REQUIRE(m_eof);	×
651
652	// Close file to trigger any exceptions due to badbit / failbit
653	try {	×
654	m_output.close();	×
655	} catch (const std::ios_base::failure&) {	×
656	std::ostringstream msg;	×
657	msg << "Error closing FASTQ file " << shell_escape(m_output.filename());	×
658
659	throw io_error(msg.str(), errno);	×
660	}	×
661	}
662
663	} // namespace adapterremoval

MikkelSchubert / adapterremoval / #110

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