Reading fasta sequence regions with faidx

Question

Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.

To do it in c++ i use the library htslib and the command faidx_fetch_seq() which is similar to read specific regions using samtools faidx.

The problem is when im reading in the reference genomes, there are of course multiple regions consisting only of NNNNN (and so on). Since i only need the actual DNA sequence i wish to read in the entire file but removing the regions with NNNN.

I have tried reading in the regions in c++ and then just creating an if-statement to check the actual bases consist of N. This works but is incredibly slow.

So my question is: do anyone know how to only read in the actual DNA sequences, removing the NNN regions? either with samtools or with the c++ funciton faidx_fetch_seq().

Answer 1

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 Ns in your FASTA file). (2) Use bedops --complement with this mask file and the output of UCSC fetchChromSizes to generate a BED file containing regions that are unmasked. (3) Repeat for other reference genomes, as needed.

You can now loop over the unmasked regions in this file, making calls to samtools faidx or faidx_fetch_seq(). These calls will return unmasked sequence.

For performance or convenience, you could read chunks of sequence into a std::deque or std::vector buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.