<p>If you can't read in the entire sequence for a given reference genome:</p>
<p>(1) Use a <a href="https://www.danielecook.com/generate-a-bedfile-of-masked-ranges-a-fasta-file/" rel="nofollow noreferrer">script</a> to generate a BED file of masked regions of low complexity (i.e. where you find <code>N</code>s in your FASTA file). (2) Use <code>bedops --complement</code> with this mask file and the output of UCSC <code>fetchChromSizes</code> to generate a BED file containing regions that are <em>unmasked</em>. (3) Repeat for other reference genomes, as needed.</p>
<p>You can now loop over the unmasked regions in this file, making calls to <code>samtools faidx</code> or <code>faidx_fetch_seq()</code>. These calls will return unmasked sequence.</p>
<p>For performance or convenience, you could read chunks of sequence into a <code>std::deque</code> or <code>std::vector</code> buffer and then read through that however many hundreds of bases, refreshing the buffer, as needed.</p>

CC BY-SA 4.0
2020-10-29T23:27:32.803Z
2020-10-29T23:27:32.803Z
Alex Reynolds
<p>Im coding in c++ and reading in different reference genomes to examine regions across the chromosomes a few hundred basepairs at a time.</p>
<p>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.</p>
<p>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.</p>
<p>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.</p>
<p>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().</p>

CC BY-SA 4.0
2020-10-29T17:03:40.017Z
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?
Bastian Schiffthaler
CC BY-SA 4.0
2020-10-29T17:35:53.723Z
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
Bastian Schiffthaler
CC BY-SA 4.0
2020-11-23T22:39:44.900Z
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.
Alex Reynolds
<samtools><phylogenetics><c++>
Reading fasta sequence regions with faidx
# Reading fasta sequence regions with faidx

(Asked by: username_0 on 2020-10-29)

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

**Tags:** samtools, phylogenetics, c++
> **Comments:**
> **username_1**: 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?  
> **username_1**: 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  
> **username_2**: 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.  

---

## 1 answers

### Answer by username_2 on 2020-10-29. Score: 1

If you can't read in the entire sequence for a given reference genome:

(1) Use a [script](https://www.danielecook.com/generate-a-bedfile-of-masked-ranges-a-fasta-file/) to generate a BED file of masked regions of low complexity (i.e. where you find `N`s 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.


---

