＃如果已安裝seqtk，請跳過以下兩行：git clone https://github.com/lh3/seqtkcd seqtk && make＃主命令行。/ seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed  

選項 -n 設置最小拉伸長度。只需按原樣使用 -p 即可。解釋它在做什麼有點複雜。

此信息存儲在序列的2bit文件表示中，因此，如果您碰巧在本地有2bit文件（或想從UCSC下載一個文件）並安裝了py2bit（您將需要版本3.0，因為我字面上是剛剛添加了對此的支持）：

  import py2bittb = py2bit.open（“ genome.2bit”）of = open（“ NNNNN.bed”， “ w”）＃一定要在tb.chroms（）。keys（）中更改chrom的輸入和輸出文件名：blocks = tb.hardMaskedBlocks（chrom），以塊為單位：of.write（“ {} \ t {} \ t {} \ n“ .format（chrom，block [0]，block [1]））。close（）tb.close（） 

的格式，如果有幫助，您可以也可以使用它來獲取所有軟掩膜區域。只需將 hardMaskedBlocks（）替換為 softMaskedBlocks（），並確保在打開文件時指定 storeMasked = True 。

使用twoBitInfo：

  $ twoBitInfo file.fa -nBed output.bed  

例如，獲取染色體上所有被N遮蓋的區域Y（還請注意，您可以使用stdout作為文件名直接寫入stdout，並使用url作為輸入，無需下載2bit文件）：

  $ twoBitInfo http：/ /hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed標準輸出| head -3chrY 0 10000chrY 44821 94821chrY 133871 222346  

您可以在此處從適用於您的操作系統的目錄中下載twoBitInfo： http://hgdownload.soe.ucsc.edu/admin / exe /

我也應該跳下去。這是我前一段時間寫的一個Perl腳本，用於以Ns的長度分割一個fasta序列：

https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit。 pl

我剛剛對其進行了修改，以吐出一個BED格式的文件，該文件顯示了標準錯誤下的N個位置。如下使用它：

  ./ fasta-nsplit.pl tig87_withN.fa 2>out.bed > out.split.fa  

輸出（BED文件） ：

  tig00000087 0 60tig00000087 900 960tig00000087 3840 3960tig00000087 14880 14940tig00000087 59520 59700tig00000087 93000 93120tig00000087 107880 107940tig00000087 109135 109140  

以下是完整的腳本，用於完整性/修復/註釋：

 ＃！/ usr / bin / perluse警告；使用嚴格；我的$ seq =“”;我的$ shortSeqID =“”;我的$ seqID =“”;我的$ keep = 0 ; my $ cumLength = 0; while（<>）{chomp; if（/ ^ >（（。+？）（。*？\ s *）？）$ /）{我的$ newID = $ 1;我的$ newShortID = $ 2; if（$ seq）{我的$ inc = 0; while（$ seq =〜s /（NNNN +）（。*）//）{我的$ nStretch = $ 1;我的$ newSeq = $ 2;如果（$ seq），printf（“ >％s。％s \ n％s \ n”，$ seqID，$ inc ++，$ seq）; $ cumLength + =長度（$ seq）; printf（STDERR“％s \ t％d \ t％d \ n”，$ shortSeqID，$ cumLength，$ cumLength + length（$ nStretch））; $ cumLength + =長度（$ nStretch）; $ seq = $ newSeq; } printf（“ >％s \ n％s \ n”，$ seqID，$ seq）如果（$ seq）; } $ seq =“”; $ shortSeqID = $ newShortID; $ seqID = $ newID; $ cumLength = 0; } else {$ seq。= $ _; }} if（$ seq）{my $ inc = 0; while（$ seq =〜s /（NNNN +）（。*）//）{我的$ nStretch = $ 1;我的$ newSeq = $ 2;如果（$ seq），printf（“ >％s。％s \ n％s \ n”，$ seqID，$ inc ++，$ seq）; $ cumLength + =長度（$ seq）; printf（STDERR“％s \ t％d \ t％d \ n”，$ shortSeqID，$ cumLength，$ cumLength + length（$ nStretch））; $ cumLength + =長度（$ nStretch）; $ seq = $ newSeq; } printf（“ >％s \ n％s \ n”，$ seqID，$ seq）if（$ seq）;}  

這是一種從基因組序列自己生成它的方法。首先，將基因組的fasta文件轉換為tbl格式（ <seq id> \ t<sequence>），然後使用perl查找所有連續的 N 序列的起始和終止位置或 n 。

  FastaToTbl hg38.fa | perl -F“ \ t” -ane'while（/ N + / ig）{for（0 .. $＃-）{print“ $ F [0] \ t $-[$ _] \ t $ + [$ _ ] \ n“}}'> hg38.n.bed  

說明

• FastaToTbl ：這是一個非常簡單的腳本，可將fasta轉換為tbl。只需將下面的行保存為 FastaToTbl 並保存在 $ PATH 中（例如〜/ bin ）中，並使其可執行即可。

   ＃！/ 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“}  

• Perl魔術。

  • -a 使 perl 的行為類似於 awk ，將每個輸入行拆分為 -F （在這種情況下為選項卡），然後將結果保存在aray @F 中。因此， $ F [0] 將是序列ID。
  • while（/ N + / ig）{：這將匹配所有連續的 N 或 n 的情況（ code> i 標誌使其不區分大小寫）。 Perl將所有匹配項的開始位置存儲在特殊數組 @-中，並將結束位置存儲在 @ + 中。
  • for（0 .. $＃-）{：遍歷從0到數組中最終索引（ $＃-）的所有數字 @-。
  • 打印“ $ F [0] \ t $-[$ _] \ t $ + [$ _] \ n” ：打印最小床格式數據。當前序列的名稱（ $ F [0] ），當前匹配項的起始位置（ $-[$ _] ）和相應的結束位置（ $ + [$ _] ）。

我只是在系統上運行了上面的代碼，並生成了 this 一個>。

• @ user172818代碼一開始就不能算N
• @terdon代碼弄錯了床的位置
• 我像下面這樣更改@terdon代碼；它對我有用

  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“; }}'