0
$\begingroup$

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().

$\endgroup$
3
  • $\begingroup$ Can you put "incredibly slow" in context? Slow is different for everyone, are we talking seconds, minutes, hours? Can you post some of your code that you wrote for this? $\endgroup$ Oct 29 '20 at 17:03
  • $\begingroup$ I just tested with a tiny program on a 400Mbp genome, reading in all the sequence while at the same time removing Ns. Took ~10.8 seconds, is that too slow? There are certainly ways to optimize this but you will be mostly bound by disk speed at this point, checking for Ns in a string of characters is not your bottleneck $\endgroup$ Oct 29 '20 at 17:35
  • $\begingroup$ Perhaps they are working with plant genomes, some of which will be unlikely to fit into the system memory of all but the most decked-out systems. $\endgroup$ Nov 23 '20 at 22:39
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.