# hvis du har seqtk installert, hopp over følgende to linjergit-klon https://github.com/lh3/seqtkcd seqtk && gjør # til hovedkommandolinjen. / seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed  

Alternativ -n angir min strekningslengde. Bare bruk -p som det er. Det er litt komplisert å forklare hva den gjør.

Denne informasjonen er lagret i 2bit-filrepresentasjonen av sekvensen, så hvis du tilfeldigvis har en 2bit-fil lokalt (eller vil laste ned en fra UCSC) og har py2bit installert (du trenger versjon 3.0, siden jeg bokstavelig talt nettopp lagt til støtte for dette):

  importer py2bittb = py2bit.open ("genom.2bit") av = åpen ("NNNNN.bed", "w") # Endre inn- og utdatafilenavn for chrom i tb.chroms (). nøkler (): blokker = tb.hardMaskedBlocks (chrom) for blokk i blokker: of.write ("{} \ t {} \ t {} \ n ".format (chrom, blokk [0], blokk [1])) av.close () tb.close ()  

Hvis det er nyttig, kan du bruk også dette for å få alle myke maskerte regioner. Bare bytt ut hardMaskedBlocks () med softMaskedBlocks () og sørg for å spesifisere storeMasked = True når du åpner filen.

Bruk twoBitInfo:

  $ twoBitInfo file.fa -nBed output.bed  

For eksempel for å få alle N-maskerte regioner på kromosom Y (merk også at du kan bruke stdout som filnavn for å skrive direkte til stdout, og bruk av url som inngang, ikke behov for å laste ned 2bit-filen):

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

Du kan laste ned twoBitInfo fra katalogen som passer til operativsystemet ditt her: http://hgdownload.soe.ucsc.edu/admin / exe /

Jeg kan like godt hoppe på også. Her er et Perl-skript som jeg skrev for en stund siden for å dele en fasta-sekvens på strekninger av Ns:

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

Jeg har nettopp modifisert den for å spytte ut en BED-formatert fil som viser N-plasseringene ved standardfeil. Bruk den som følger:

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

Output (BED file) :

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

Her er full scripting / for fullstendighet :

  #! / usr / bin / perluse advarsler; bruk streng; min $ seq = ""; min $ shortSeqID = ""; min $ seqID = ""; min $ keep = 0 ; min $ cumLength = 0; mens (<>) {chomp; hvis (/ ^ > ((. +?) (. *? \ s *)?) $ /) {min $ newID = $ 1; min $ newShortID = $ 2; hvis ($ seq) {my $ inc = 0; mens ($ seq = ~ s / (NNNN +) (. *) //) {min $ nStretch = $ 1; min $ newSeq = $ 2; printf (">% s.% s \ n% s \ n", $ seqID, $ inc ++, $ seq) if ($ seq); $ cumLength + = lengde ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + lengde ($ nStretch)); $ cumLength + = lengde ($ nStretch); $ seq = $ newSeq; } printf (">% s \ n% s \ n", $ seqID, $ seq) if ($ seq); } $ seq = ""; $ shortSeqID = $ newShortID; $ seqID = $ newID; $ cumLength = 0; } annet {$ seq. = $ _; }} hvis ($ seq) {my $ inc = 0; mens ($ seq = ~ s / (NNNN +) (. *) //) {min $ nStretch = $ 1; min $ newSeq = $ 2; printf (">% s.% s \ n% s \ n", $ seqID, $ inc ++, $ seq) if ($ seq); $ cumLength + = lengde ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + lengde ($ nStretch)); $ cumLength + = lengde ($ nStretch); $ seq = $ newSeq; } printf (">% s \ n% s \ n", $ seqID, $ seq) if ($ seq);}  

Her er en måte å generere det selv fra genomssekvensen. Konverter først fasta-filen til genomet til tbl-format ( <seq id> \ t<sequence>) , og bruk deretter perl for å finne start- og sluttposisjonen til alle strekninger av påfølgende av N eller n .

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

Forklaring

• FastaToTbl : dette er en veldig enkelt skript som konverterer fasta til tbl. Bare lagre linjene nedenfor et sted i $ PATH (f.eks. ~ / bin ) som FastaToTbl og gjør den kjørbar.

    #! / usr / bin / awk -f {if (substr ($ 1,1,1) == ">") if (NR>1) printf "\ n% s \ t" , substr ($ 0,2, lengde ($ 0) -1) annet printf "% s \ t", substr ($ 0,2, lengde ($ 0) -1) annet printf "% s", $ 0} END {printf "\ n "}  

• Perl-magien.

  • -a gjør at perl oppfører seg som awk , og deler hver inngangslinje på verdien gitt av -F (en fane i dette tilfellet) og lagre resultatet i den aray @F . Så, $ F [0] vil være sekvens-ID.
  • while (/ N + / ig) {: dette vil matche alle tilfeller av påfølgende N eller n ( i flagg gjør det ikke skiftende på store og små bokstaver). Perl vil lagre startposisjonene til alle treff i den spesielle matrisen @ - og sluttposisjonene i @ + .
  • for (0 .. $ # -) {: gjenta over alle tall fra 0 til den endelige indeksen ( $ # - ) i matrisen @ - .
  • skriv ut "$ F [0] \ t $ - [$ _] \ t $ + [$ _] \ n" : skriv ut minimale sengeformatdata. Navnet på den gjeldende sekvensen ( $ F [0] ), startposisjonen til den gjeldende kampen ( $ - [$ _] ) og den tilsvarende sluttposisjonen ( $ + [$ _] ).

Jeg kjørte nettopp over på systemet mitt, og det genererte dette.

• @ user172818-koden kan ikke telle N fra begynnelsen
• @terdon-kode får feil sengeposisjon
• jeg endrer @terdon-kode som nedenfor; og det fungerer for meg

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