Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-10-29T23:27:32.803Z
Answers.LastActivityDate: 2020-10-29T23:27:32.803Z
Answers.UserDisplayName: Alex Reynolds
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-10-29T17:03:40.017Z
Comments.Text: 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?
Comments.UserDisplayName: Bastian Schiffthaler
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-10-29T17:35:53.723Z
Comments.Text: 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
Comments.UserDisplayName: Bastian Schiffthaler
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-11-23T22:39:44.900Z
Comments.Text: 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.
Comments.UserDisplayName: Alex Reynolds
Tags: <samtools><phylogenetics><c++>
Title: Reading fasta sequence regions with faidx
ThreadText: # 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.
  
  
  ---
  
  