<p>This information is stored in the 2bit file representation of the sequence, so if you happen to have a 2bit file locally (or want to download one from UCSC) and have py2bit installed (you'll need version 3.0, since I literally just added support for this):</p>

<pre class="lang-py prettyprint-override"><code>import py2bit
tb = py2bit.open("genome.2bit")
of = open("NNNNN.bed", "w")  # Do change the input and output file names
for chrom in tb.chroms().keys():
    blocks = tb.hardMaskedBlocks(chrom)
    for block in blocks:
        of.write("{}\t{}\t{}\n".format(chrom, block[0], block[1]))
of.close()
tb.close()
</code></pre>

<p>In case it's helpful, you can also use this to get all soft-masked regions. Just replace <code>hardMaskedBlocks()</code> with <code>softMaskedBlocks()</code> and ensure to specify <code>storeMasked=True</code> when you open the file.</p>

CC BY-SA 3.0
2018-02-13T12:16:07.913Z
2018-02-13T12:48:30.077Z
2018-02-13T12:48:30.077Z
Devon Ryan
<p>Use twoBitInfo:</p>

<pre><code>$ twoBitInfo file.fa -nBed output.bed
</code></pre>

<p>For example, to get all the N-masked regions on chromosome Y (also note you can use stdout as a filename to write directly to stdout, and use of the url as the input, no need to download the 2bit file):</p>

<pre><code>$ twoBitInfo http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed stdout | head -3
chrY    0   10000
chrY    44821   94821
chrY    133871  222346
</code></pre>

<p>You can download twoBitInfo from the directory appropriate to your operating system here:
<a href="http://hgdownload.soe.ucsc.edu/admin/exe/" rel="nofollow noreferrer">http://hgdownload.soe.ucsc.edu/admin/exe/</a></p>

CC BY-SA 3.0
2018-02-14T19:58:31.830Z
2018-02-14T19:58:31.830Z
Christopher Lee
<p>I might as well jump on as well. Here's a Perl script that I wrote a while ago to split a fasta sequence at stretches of Ns:</p>

<p><a href="https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit.pl" rel="nofollow noreferrer">https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit.pl</a></p>

<p>I've just modified it to spit out a BED formatted file showing the N locations on standard error. Use it as follows:</p>

<pre><code>./fasta-nsplit.pl tig87_withN.fa 2&gt;out.bed &gt; out.split.fa
</code></pre>

<p>Output (BED file):</p>

<pre><code>tig00000087 0   60
tig00000087 900 960
tig00000087 3840    3960
tig00000087 14880   14940
tig00000087 59520   59700
tig00000087 93000   93120
tig00000087 107880  107940
tig00000087 109135  109140
</code></pre>

<p>Here's the full script, for completeness / fixing / comment:</p>

<pre><code>#!/usr/bin/perl
use warnings;
use strict;

my $seq = "";
my $shortSeqID = "";
my $seqID = "";
my $keep = 0;
my $cumLength = 0;
while(&lt;&gt;){
  chomp;
  if(/^&gt;((.+?)( .*?\s*)?)$/){
    my $newID = $1;
    my $newShortID = $2;
    if($seq){
      my $inc = 0;
      while($seq =~ s/(NNNN+)(.*)//){
 my $nStretch = $1;
        my $newSeq = $2;
        printf("&gt;%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
 $cumLength += length($seq);
 printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength, 
        $cumLength + length($nStretch));
 $cumLength += length($nStretch);
        $seq = $newSeq;
      }
      printf("&gt;%s\n%s\n", $seqID, $seq) if ($seq);
    }
    $seq = "";
    $shortSeqID = $newShortID;
    $seqID = $newID;
    $cumLength = 0;
  } else {
    $seq .= $_;
  }
}
if($seq){
  my $inc = 0;
  while($seq =~ s/(NNNN+)(.*)//){
    my $nStretch = $1;
    my $newSeq = $2;
    printf("&gt;%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
    $cumLength += length($seq);
    printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength, 
       $cumLength + length($nStretch));
    $cumLength += length($nStretch);
    $seq = $newSeq;
  }
  printf("&gt;%s\n%s\n", $seqID, $seq) if ($seq);
}
</code></pre>

CC BY-SA 3.0
2018-02-13T19:17:43.693Z
2018-02-14T03:51:29.297Z
2018-02-14T03:51:29.297Z
gringer
<ul>
<li>@user172818 code can't count N at fa begining</li>
<li>@terdon code get wrong bed position</li>
<li>i change @terdon code like below; and it works for me</li>
</ul>

<pre><code>awk '{
    if (substr($1,1,1)=="&gt;")
        if (NR&gt;1)
            printf "\n%s,", substr($0,2,length($0)-1)
        else
            printf "%s,", substr($0,2,length($0)-1)
    else
        printf "%s", $0
}
END{printf "\n"}' test.fa | \
perl -F"," -ane 'while($F[1] =~ /N+/ig){
                        for(0..$#-){                    
                            print "$F[0]\t$-[$_]\t$+[$_]\n";
                        }
                  }'
</code></pre>

CC BY-SA 4.0
2019-11-01T02:45:43.000Z
2019-11-01T07:46:25.240Z
2019-11-01T07:46:25.240Z
atongsa
<p>Here's a way to generate it yourself from the genome sequence. First, convert the fasta file of the genome to tbl format (<code>&lt;seq id&gt;\t&lt;sequence&gt;)</code>, then use perl to find the start and end positions of all stretches of consecutive of <code>N</code> or <code>n</code>. </p>

<pre><code>FastaToTbl hg38.fa | 
    perl -F"\t" -ane 'while(/N+/ig){
                            for(0..$#-){
                                print "$F[0]\t$-[$_]\t$+[$_]\n"
                            }
                      }' &gt; hg38.n.bed
</code></pre>

<h3>Explanation</h3>

<ul>
<li><p><code>FastaToTbl</code> : this is a very simple script that converts fasta to tbl. Just save the lines below somewhere in your <code>$PATH</code> (e.g. <code>~/bin</code>) as <code>FastaToTbl</code> and make it executable. </p>

<pre><code>#!/usr/bin/awk -f
{
    if (substr($1,1,1)=="&gt;")
            if (NR&gt;1)
                printf "\n%s\t", substr($0,2,length($0)-1)
            else 
                printf "%s\t", substr($0,2,length($0)-1)
        else 
            printf "%s", $0
}
END{printf "\n"}
</code></pre></li>
<li><p>The Perl magic. </p>

<ul>
<li>The <code>-a</code> makes <code>perl</code> behave like <code>awk</code>, splitting each input line on the value given by <code>-F</code> (a tab in this case) and saving the result in the aray <code>@F</code>. So, <code>$F[0]</code> will be the sequence id. </li>
<li><code>while(/N+/ig){</code> : this will match all cases of consecutive <code>N</code> or <code>n</code> (the <code>i</code> flag makes it case-insensitive). Perl will store the start positions of all matches in the special array <code>@-</code> and the end positions in <code>@+</code>. </li>
<li><code>for(0..$#-){</code> : iterate over all numbers from 0 to the final index (<code>$#-</code>) in the array <code>@-</code>.</li>
<li><code>print "$F[0]\t$-[$_]\t$+[$_]\n"</code> : print minimal bed format data. The name of the current sequence (<code>$F[0]</code>), the start position of the current match (<code>$-[$_]</code>) and the corresponding end position (<code>$+[$_]</code>).</li>
</ul></li>
</ul>

<p>I just ran the above on my system and it generated <a href="https://pastebin.com/UbakGp9b" rel="nofollow noreferrer">this</a>. </p>

CC BY-SA 3.0
2018-02-13T10:53:30.140Z
2018-02-13T10:53:30.140Z
terdon
<pre class="lang-sh prettyprint-override"><code># if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk &amp;&amp; make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa &gt; hg38-N.bed
</code></pre>

<p>Option <code>-n</code> sets the min stretch length. Just use <code>-p</code> as is. It is a bit complicated to explain what it is doing.</p>

CC BY-SA 3.0
2018-02-13T14:25:04.800Z
2018-02-13T14:25:04.800Z
user172818
<p>How do I obtain a bed file with the list of N-regions of GRCh38 reference? This is, the regions where the sequence is a stretch of Ns.</p>

<genome><bed>
bed file with N regions of GRCh38 reference?
# bed file with N regions of GRCh38 reference?

(Asked by: username_0 on 2018-02-13)

How do I obtain a bed file with the list of N-regions of GRCh38 reference? This is, the regions where the sequence is a stretch of Ns.

**Tags:** genome, bed

---

## 6 answers

### Answer by username_6 on 2018-02-13. Score: 5 | Accepted answer

# if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk && make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed

Option -n sets the min stretch length. Just use -p as is. It is a bit complicated to explain what it is doing.


---

### Answer by username_1 on 2018-02-13. Score: 4

This information is stored in the 2bit file representation of the sequence, so if you happen to have a 2bit file locally (or want to download one from UCSC) and have py2bit installed (you'll need version 3.0, since I literally just added support for this):
import py2bit
tb = py2bit.open("genome.2bit")
of = open("NNNNN.bed", "w")  # Do change the input and output file names
for chrom in tb.chroms().keys():
    blocks = tb.hardMaskedBlocks(chrom)
    for block in blocks:
        of.write("{}\t{}\t{}\n".format(chrom, block[0], block[1]))
of.close()
tb.close()

In case it's helpful, you can also use this to get all soft-masked regions. Just replace hardMaskedBlocks() with softMaskedBlocks() and ensure to specify storeMasked=True when you open the file.


---

### Answer by username_2 on 2018-02-14. Score: 3

Use twoBitInfo:
$ twoBitInfo file.fa -nBed output.bed

For example, to get all the N-masked regions on chromosome Y (also note you can use stdout as a filename to write directly to stdout, and use of the url as the input, no need to download the 2bit file):
$ twoBitInfo http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed stdout | head -3
chrY    0   10000
chrY    44821   94821
chrY    133871  222346

You can download twoBitInfo from the directory appropriate to your operating system here:
http://hgdownload.soe.ucsc.edu/admin/exe/


---

### Answer by username_3 on 2018-02-13. Score: 2

I might as well jump on as well. Here's a Perl script that I wrote a while ago to split a fasta sequence at stretches of Ns:
https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit.pl
I've just modified it to spit out a BED formatted file showing the N locations on standard error. Use it as follows:
./fasta-nsplit.pl tig87_withN.fa 2>out.bed > out.split.fa

Output (BED file):
tig00000087 0   60
tig00000087 900 960
tig00000087 3840    3960
tig00000087 14880   14940
tig00000087 59520   59700
tig00000087 93000   93120
tig00000087 107880  107940
tig00000087 109135  109140

Here's the full script, for completeness / fixing / comment:
#!/usr/bin/perl
use warnings;
use strict;

my $seq = "";
my $shortSeqID = "";
my $seqID = "";
my $keep = 0;
my $cumLength = 0;
while(<>){
  chomp;
  if(/^>((.+?)( .*?\s*)?)$/){
    my $newID = $1;
    my $newShortID = $2;
    if($seq){
      my $inc = 0;
      while($seq =~ s/(NNNN+)(.*)//){
 my $nStretch = $1;
        my $newSeq = $2;
        printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
 $cumLength += length($seq);
 printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength, 
        $cumLength + length($nStretch));
 $cumLength += length($nStretch);
        $seq = $newSeq;
      }
      printf(">%s\n%s\n", $seqID, $seq) if ($seq);
    }
    $seq = "";
    $shortSeqID = $newShortID;
    $seqID = $newID;
    $cumLength = 0;
  } else {
    $seq .= $_;
  }
}
if($seq){
  my $inc = 0;
  while($seq =~ s/(NNNN+)(.*)//){
    my $nStretch = $1;
    my $newSeq = $2;
    printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
    $cumLength += length($seq);
    printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength, 
       $cumLength + length($nStretch));
    $cumLength += length($nStretch);
    $seq = $newSeq;
  }
  printf(">%s\n%s\n", $seqID, $seq) if ($seq);
}


---

### Answer by username_4 on 2019-11-01. Score: 1

@user172818 code can't count N at fa begining
@terdon code get wrong bed position
i change @terdon code like below; and it works for me

awk '{
    if (substr($1,1,1)==">")
        if (NR>1)
            printf "\n%s,", substr($0,2,length($0)-1)
        else
            printf "%s,", substr($0,2,length($0)-1)
    else
        printf "%s", $0
}
END{printf "\n"}' test.fa | \
perl -F"," -ane 'while($F[1] =~ /N+/ig){
                        for(0..$#-){                    
                            print "$F[0]\t$-[$_]\t$+[$_]\n";
                        }
                  }'


---

### Answer by username_5 on 2018-02-13. Score: 1

Here's a way to generate it yourself from the genome sequence. First, convert the fasta file of the genome to tbl format (\t), then use perl to find the start and end positions of all stretches of consecutive of N or n. 
FastaToTbl hg38.fa | 
    perl -F"\t" -ane 'while(/N+/ig){
                            for(0..$#-){
                                print "$F[0]\t$-[$_]\t$+[$_]\n"
                            }
                      }' > hg38.n.bed

Explanation

FastaToTbl : this is a very simple script that converts fasta to tbl. Just save the lines below somewhere in your $PATH (e.g. ~/bin) as FastaToTbl and make it executable. 
#!/usr/bin/awk -f
{
    if (substr($1,1,1)==">")
            if (NR>1)
                printf "\n%s\t", substr($0,2,length($0)-1)
            else 
                printf "%s\t", substr($0,2,length($0)-1)
        else 
            printf "%s", $0
}
END{printf "\n"}

The Perl magic. 

The -a makes perl behave like awk, splitting each input line on the value given by -F (a tab in this case) and saving the result in the aray @F. So, $F[0] will be the sequence id. 
while(/N+/ig){ : this will match all cases of consecutive N or n (the i flag makes it case-insensitive). Perl will store the start positions of all matches in the special array @- and the end positions in @+. 
for(0..$#-){ : iterate over all numbers from 0 to the final index ($#-) in the array @-.
print "$F[0]\t$-[$_]\t$+[$_]\n" : print minimal bed format data. The name of the current sequence ($F[0]), the start position of the current match ($-[$_]) and the corresponding end position ($+[$_]).


I just ran the above on my system and it generated this.


---

