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}$...
- 6,485
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 ...
- 3,125
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 ...
- 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 ...
- 6,485
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 ...
- 6,485
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 ...
- 13.8k
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 ...
- 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 ...
- 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 ...
- 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 ...
- 956
4
votes
Accepted
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 ...
- 189
3
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 ...
- 19.6k
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....
- 4,845
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 ...
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 ...
- 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 ...
- 2,739
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-...
- 6,485
2
votes
Accepted
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...
- 11
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/...
- 13.8k
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,...,...
- 810
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
...
- 1,501
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 ...
- 3,125
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 ...
Only top scored, non community-wiki answers of a minimum length are eligible