Skip to main content
7 votes
Accepted

How to write a hash function for canonical kmers?

You can convert any string hash function to a "canonical" DNA string hash function. Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$...
user172818's user avatar
  • 6,535
7 votes
Accepted

How to calculate the memory usage of storing kmers in RAM

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table. Note: I wrote a C++ port ...
Alex Reynolds's user avatar
7 votes
Accepted

Multithread fastq processing with kseq.h in C++11?

You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the ...
nomad's user avatar
  • 472
6 votes
Accepted

What is the most compact data structure for canonical k-mers with the fastest lookup time?

The question and the accepted answer are not about k-mer data structure at all, which I will explain in detail below. I will first answer the actual question OP intends to ask. The simplest way to ...
user172818's user avatar
  • 6,535
6 votes

How to calculate the memory usage of storing kmers in RAM

I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM? Hash table based k-mer counters only keep k-mers seen in data. For $16<k\le32$, you need 64-bit ...
user172818's user avatar
  • 6,535
5 votes

How to write a hash function for canonical kmers?

However you want. One way to do this is to generate both the forward and reverse complement kmers, then choose the lexicographically-least kmer for the storage key. To delve further into this ...
gringer's user avatar
  • 14.3k
5 votes

How to calculate the memory usage of storing kmers in RAM

If you simply need a set of present/absent k-mers and your kmer length is short, then convert the sequence length to an integer (a 13-mer would require 26 bits or 4 bytes) and use that as an index ...
Devon Ryan's user avatar
  • 19.6k
4 votes
Accepted

How to build a BAM header file with htslib in C++?

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there). With BAM files one doesn't ...
Devon Ryan's user avatar
  • 19.6k
4 votes

Read sorted and indexed BAM files faster in C++?

In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest. In the first method I use an ...
Nisba's user avatar
  • 141
4 votes
Accepted

Use Htslib to extract specific region of vcf file

For iterating over reads in a region in a single BCF file, you want to use the iterator interface: bcf_itr_queryi() or ...
John Marshall's user avatar
4 votes
Accepted

How do I access sequence and MAPQ from bam using HTSlib in C++?

The MAPQ is stored in b->core.qual, given a bam1_t *b. You can access the sequence with ...
Devon Ryan's user avatar
  • 19.6k
4 votes
Accepted

How to use C++ htslib to read VCF contig name and size?

...
user172818's user avatar
  • 6,535
3 votes

Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'

Htslib is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM: ...
3 votes
Accepted

Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'

The helpdesk for my institution's supercomputers showed me the following flags: g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread This fixed the error. I'm now ...
annabelperry's user avatar
3 votes

While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError'

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message....
Konrad Rudolph's user avatar
3 votes

Read sorted and indexed BAM files faster in C++?

If you need to iterate through the entire BAM file, decompression is usually the main bottleneck. Three things you can do: Install libdeflate (https://github.com/ebiggers/libdeflate ) and then the ...
Christopher Chang's user avatar
2 votes

Read sorted and indexed BAM files faster in C++?

In both deepTools and MethylDackel we implement multithreading by having individual threads process separate genomic regions. At some point you hit an IO limit, but unless your processing is trivial ...
Devon Ryan's user avatar
  • 19.6k
2 votes
Accepted

Non-random access on a fastq file

You don't have a RAM problem because you only have about: 10-20 reads that I want to access from each file. Please go for (2). If you are concerned about RAM usage, you may selectively write some ...
SmallChess's user avatar
  • 2,699
2 votes
Accepted

How is consensus alignment for OLC assembly usually implemented?

For short reads, mainstream assemblers don't have a consensus step. The contig path directly spells the consensus sequence. For long reads, the most popular way to derive a consensus is via partial-...
user172818's user avatar
  • 6,535
2 votes
Accepted

Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?

...
user172818's user avatar
  • 6,535
1 vote

Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'

I think SEQAN_HAS_ZLIB flag not works correctly as it supposed to be. So you can check the seqan library to see if SEQAN_HAS_ZLIB...
xin king's user avatar
1 vote

Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'

I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the SeqAn BAM tutorial, it looks like you need to use the BamFileIn object for reading BAM/...
gringer's user avatar
  • 14.3k
1 vote

Construct the Overlap Graph of a Collection of k-mers

You will want to create a self loop for a $k$-mer $s$ whenever $s_1,...,s_{k-1} = s_2,...,s_k$. This will only happen when all of the characters of the $k$-mer are identical, as $s_1=s_2, s_2=s_3,...,...
Throckmorton's user avatar
1 vote

How to print alternative allele from VCF with htslib c++

see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104 ...
Pierre's user avatar
  • 1,536
1 vote

Reading fasta sequence regions with faidx

If you can't read in the entire sequence for a given reference genome: (1) Use a script to generate a BED file of masked regions of low complexity (i.e. where you find ...
Alex Reynolds's user avatar
1 vote

c++ program to find length string an char*

I think you're confusing C++ and C strings, which are very different. When you have a char *, that's a C-style "string". It's really just a bunch of bytes with a ...
Jose Fernando Lopez Fernandez's user avatar

Only top scored, non community-wiki answers of a minimum length are eligible