Skip to content

Allow for correct counting of FASTQ records which have quality string beginning with '@' #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ahstram wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Allow for correct counting of FASTQ records which have quality string beginning with '@' #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1ebd9e7
Remove trailing whitespaces
Feb 16, 2019
e0341f9
Allow for correct counting of FASTQ records which have quality string…
Feb 16, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 35 additions & 31 deletions bioio.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
/* bioio.hpp -- FASTA/Q I/O

Copyright (C) 2017 University of Oxford.

Author: Daniel Cooke <dcooke@well.ox.ac.uk>

Use of this source code is governed by the MIT license that can be found in the LICENSE file. */

#ifndef __bioio__bioio__
Expand Down Expand Up @@ -48,19 +48,19 @@ struct FastaContigIndex
std::size_t offset;
std::size_t line_length;
std::size_t line_byte_length;

FastaContigIndex() = default;

template <typename T>
explicit FastaContigIndex(T&& contig_name, std::size_t length, std::size_t offset,
std::size_t line_length, std::size_t line_byte_length)
: contig_name {std::forward<T>(contig_name)}
, offset {offset}
, length {length}
, line_length {line_length}
, line_byte_length {line_byte_length}
, line_byte_length {line_byte_length}
{}

template <typename T>
explicit FastaContigIndex(const T& fasta_index_line)
{
Expand All @@ -80,7 +80,7 @@ struct FastaRecord
{
StringType name;
SequenceType sequence;

FastaRecord() = delete;
template <typename StringType_, typename SequenceType_>
explicit FastaRecord(StringType_&& name, SequenceType_&& sequence)
Expand All @@ -96,13 +96,13 @@ struct FastqRecord
StringType name;
SequenceType1 seq;
SequenceType2 qual;

FastqRecord() = delete;
template <typename StringType_, typename SequenceType1_, typename SequenceType2_>
explicit FastqRecord(StringType_&& name, SequenceType1_&& seq, SequenceType2_&& qual)
:name {std::forward<StringType_>(name)}
, seq {std::forward<SequenceType1_>(seq)}
, qual {std::forward<SequenceType2_>(qual)}
, qual {std::forward<SequenceType2_>(qual)}
{}
};

Expand All @@ -125,23 +125,27 @@ using FastqReads = std::pair<ReadIds<StringType>, FastqMap<StringType, SequenceT

namespace detail
{
static constexpr char fasta_delim {'>'};
static constexpr char fastq_delim {'@'};

// Counts the occurrences of record_delim that follow a newline '\n'
inline size_t count_records(std::istream& is, const char record_delim)
{
const auto current_position = is.tellg();
size_t result {};
while (is) {
if (is.peek() == record_delim) ++result;
is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
if (record_delim == fastq_delim) {
result = std::count(std::istreambuf_iterator<char>(is), std::istreambuf_iterator<char>(), '\n')/4;
} else {
while (is) {
if (is.peek() == record_delim) ++result;
is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}
}
is.clear();
is.seekg(current_position, std::ios::beg);
return result;
}

static constexpr char fasta_delim {'>'};
static constexpr char fastq_delim {'@'};


template <typename StringType = std::string, typename SequenceType = std::string>
::bioio::FastaRecord<StringType, SequenceType>
read_fasta_record(std::istream& fasta)
Expand All @@ -150,7 +154,7 @@ namespace detail
std::getline(fasta, name); // name is always a single line
SequenceType line;
std::getline(fasta, line);

// The FASTA format is not as simple as FASTQ - the sequence
// may be broken into multiple lines. We assume each line is the same size.
if (!fasta.good() || fasta.peek() == fasta_delim) {
Expand Down Expand Up @@ -178,7 +182,7 @@ namespace detail
return ::bioio::FastaRecord<StringType, SequenceType> {std::move(name), std::move(seq)};
}
}

template <typename StringType = std::string, typename SequenceType1 = std::string,
typename SequenceType2 = std::string>
::bioio::FastqRecord<StringType, SequenceType1, SequenceType2>
Expand All @@ -187,14 +191,14 @@ namespace detail
StringType name;
SequenceType1 seq;
SequenceType2 quals;

// Unlike FASTA, FASTQ always use one line per field.
std::getline(fastq, name);
std::getline(fastq, seq);
fastq.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
std::getline(fastq, quals);
return ::bioio::FastqRecord<StringType, SequenceType1, SequenceType2>

return ::bioio::FastqRecord<StringType, SequenceType1, SequenceType2>
{std::move(name), std::move(seq), std::move(quals)};
}
} // namespace detail
Expand All @@ -205,7 +209,7 @@ namespace detail

inline std::size_t count_contigs_in_fasta_index(std::istream& fasta_index)
{
return std::count(std::istreambuf_iterator<char>(fasta_index),
return std::count(std::istreambuf_iterator<char>(fasta_index),
std::istreambuf_iterator<char>(), '\n');
}

Expand Down Expand Up @@ -296,7 +300,7 @@ inline std::size_t remaining_line_length(const ::bioio::FastaContigIndex& index,
} // namespace detail

template <typename SequenceType = std::string>
SequenceType
SequenceType
read_fasta_contig(std::istream& fasta, const FastaContigIndex& index,
const std::size_t begin, std::size_t length)
{
Expand Down Expand Up @@ -613,7 +617,7 @@ template <typename StringType = std::string, typename SequenceType = std::string
FastaReads<StringType, SequenceType>
read_fasta_map(const std::string& path, const ReadIds<StringType>& names)
{
return read_fasta_map(path, names,
return read_fasta_map(path, names,
[] (StringType&& name) { return name; });
}

Expand Down Expand Up @@ -715,7 +719,7 @@ std::vector<FastqRecord<StringType, SequenceType1, SequenceType2>>
read_fastq(const std::string& path)
{
return read_fastq<StringType, SequenceType1, SequenceType2>(path, [] (const StringType& name) { return name; });
}
}

template <typename StringType = std::string, typename SequenceType1 = std::string,
typename SequenceType2 = std::string>
Expand All @@ -730,8 +734,8 @@ template <typename StringType = std::string, typename SequenceType1 = std::strin
std::vector<FastqRecord<StringType, SequenceType1, SequenceType2>>
read_fastq(std::istream& fastq, const std::size_t num_records)
{
return read_fastq<StringType, SequenceType1, SequenceType2>(fastq, num_records, [] (const StringType& name) {
return name;
return read_fastq<StringType, SequenceType1, SequenceType2>(fastq, num_records, [] (const StringType& name) {
return name;
});
}

Expand All @@ -741,8 +745,8 @@ std::vector<FastqRecord<StringType, SequenceType1, SequenceType2>>
read_fastq(const std::string& path, const std::size_t num_records)
{
std::ifstream fastq {path, std::ios::binary};
return read_fastq<StringType, SequenceType1, SequenceType2>(fastq, num_records, [] (const StringType& name) {
return name;
return read_fastq<StringType, SequenceType1, SequenceType2>(fastq, num_records, [] (const StringType& name) {
return name;
});
}

Expand All @@ -763,7 +767,7 @@ read_fastq_map(const std::string& path)
}
return {names, data};
}

} // namespace bioio

#endif
56 changes: 28 additions & 28 deletions fasta.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
/* fasta.cpp

Copyright (C) 2017 University of Oxford.

Author: Daniel Cooke <dcooke@well.ox.ac.uk>

Use of this source code is governed by the MIT license that can be found in the LICENSE file */

#include <iostream>
Expand All @@ -22,51 +22,51 @@ using GenomicRegion = std::tuple<std::string, size_t, size_t>;
GenomicRegion parse_region(std::string region, const bioio::FastaIndex& index)
{
region.erase(std::remove(region.begin(), region.end(), ','), region.end());

const static std::regex re {"([^:]+)(?::(\\d+)(-)?(\\d*))?"};

std::smatch match;

if (std::regex_match(region, match, re) && match.size() == 5) {
auto contig_name = match.str(1);

if (index.count(contig_name) == 0) {
throw std::runtime_error {"contig " + contig_name + " not found"};
}

const auto contig_size = index.at(contig_name).length;

size_t begin {0}, end {0};

if (match.str(2).empty()) {
end = contig_size;
} else {
begin = static_cast<size_t>(std::stoull(match.str(2)));

if (match.str(3).empty()) {
end = begin + 1;
} else if (match.str(4).empty()) {
end = contig_size;
} else {
end = static_cast<size_t>(std::stoull(match.str(4)));
}

if (begin > contig_size) {
throw std::runtime_error {"region " + region + " is larger than contig " + contig_name + ":0-" + std::to_string(contig_size)};
}

if (begin > end) {
throw std::runtime_error {"begin position is past end position in region " + region};
}

if (end > contig_size) {
end = contig_size;
}
}

return GenomicRegion {std::move(contig_name), begin, end};
}

throw std::runtime_error {"could not parse region " + region};
}

Expand Down Expand Up @@ -101,50 +101,50 @@ int main(int argc, char **argv)
print_usage();
return 1;
}

auto index_path = get_cmd_option(argv, argv + argc, "-i");

std::string fasta_path {argv[argc - 2]};

std::ifstream fasta {fasta_path, std::ios::binary};

if (!fasta) {
std::cerr << "Error: could not open fasta " << fasta_path << std::endl;
print_usage();
return 1;
}

if (index_path.empty() && index_path.find(".") != std::string::npos) {
index_path = fasta_path;
index_path.replace(index_path.begin() + index_path.find_last_of("."), index_path.end(), ".fai");
}

std::ifstream index_file {index_path, std::ios::binary};

if (!index_file) {
index_file.open(fasta_path + ".fai");
if (!index_file) {
std::cerr << "Error: could not open index file, use samtools faidx <fasta> to make index file" << std::endl;
return 1;
}
}

try {
const auto index = bioio::read_fasta_index(index_file);

const auto region = parse_region(argv[argc - 1], index);
const auto& contig = std::get<0>(region);
const auto begin = std::get<1>(region);
const auto length = std::get<2>(region) - begin;

if (cmd_option_exists(argv, argv + argc, "-s")) {
std::cout << length << std::endl;
} else {
const auto sequence = bioio::read_fasta_contig(fasta, index.at(contig), begin, length);

std::cout << sequence << std::endl;
}

return 0;
} catch (std::runtime_error& e) {
std::cerr << "Error: " << e.what() << std::endl;
Expand Down
Loading