# ha a seqtk telepítve van, hagyja ki a következő két linegit klónt: https://github.com/lh3/seqtkcd seqtk && tegye a fő parancssorba. / seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed  

A -n opció beállítja a minimális nyújtási hosszat. Csak használja a -p -t úgy, ahogy van. Kicsit bonyolult elmagyarázni, mit csinál.

Ez az információ a szekvencia 2 bites fájlábrázolásában van tárolva, tehát ha előfordul, hogy helyben van egy 2 bites fájl (vagy le akarja tölteni az UCSC-ből), és telepítve van a py2bit (a 3.0 verzióra lesz szükség, mivel én szó szerint csak hozzáadta ennek a támogatását):

  import py2bittb = py2bit.open ("genome.2bit") of = open ("NNNNN.bed", "w") # Változtassa meg a bemeneti és kimeneti fájlneveket a chrom fájlhoz a tb.chroms (). kulcsokban (): blocks = tb.hardMaskedBlocks (chrom) a következő blokkok blokkjához: t {} \ n ".formátum (chrom, [0], [1] blokk)) a.close () tb.close ()  

Ha hasznos, használja ezt az összes lágyruhás régió megszerzéséhez is. Csak cserélje le a hardMaskedBlocks () -ot a softMaskedBlocks () -ra és győződjön meg róla, hogy a fájl megnyitásakor megadja-e a storeMasked = True -ot.

A twoBitInfo használata:

  $ twoBitInfo file.fa -nBed output.bed  

Például az összes N-maszkolt régió kromoszómára jutásához Y (azt is vegye figyelembe, hogy az stdout fájlnévként közvetlenül írhatja az stdout-ra, és az URL-t használhatja bemenetként, nincs szükség a 2 bites fájl letöltésére):

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

A twoBitInfo-t letöltheti az operációs rendszernek megfelelő könyvtárból: http://hgdownload.soe.ucsc.edu/admin / exe /

Lehet, hogy én is továbbugrom. Itt van egy Perl-szkript, amelyet egy ideje írtam egy Fasta-szekvencia felosztására az N-ek szakaszain:

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

Most módosítottam, hogy kiköpjön egy BED formátumú fájlt, amely az N helyet mutatja standard hibán. Használja a következőképpen:

  ./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 107880 107940tig00000087 109135 101340 :  
  #! / usr / bin / perluse figyelmeztetések; használja szigorúan; 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; míg ($ seq = ~ s / (NNNN +) (. *) //) {my $ nStretch = $ 1; my $ newSeq = $ 2; printf (">% s.% s \ n% s \ n", $ seqID, $ inc ++, $ seq) if ($ seq); $ cumLength + = hossz ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + length ($ nStretch)); $ cumLength + = hossz ($ 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; míg ($ seq = ~ s / (NNNN +) (. *) //) {my $ nStretch = $ 1; my $ newSeq = $ 2; printf (">% s.% s \ n% s \ n", $ seqID, $ inc ++, $ seq) if ($ seq); $ cumLength + = hossz ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + length ($ nStretch)); $ cumLength + = hossz ($ nStretch); $ seq = $ newSeq; } printf (">% s \ n% s \ n", $ seqID, $ seq) if ($ seq);}  

Itt van egy módszer arra, hogy saját maga állítsa elő a genomszekvenciából. Először konvertálja a genom fasta fájlját tbl formátumra ( <seq id> \ t<sequence>) , majd a perl használatával keresse meg a N összes egymást követő szakaszának kezdő és befejező pozícióját vagy n .

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

Magyarázat

• FastaToTbl : ez egy nagyon egyszerű szkript, amely a fastát tbl-vé alakítja. Csak mentse az alábbi sorokat valahova a $ PATH -ba (pl. ~ / bin ) FastaToTbl néven, és futtathatóvá tegye.

    #! / 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 "}  

• A Perl varázslat.

  • A -a hatására a perl úgy viselkedik, mint az awk , az egyes bemeneti sorokat elosztva az -F (ebben az esetben egy fül), és az eredményt elmenti az aray @F fájlba. Tehát a $ F [0] lesz a sorozat azonosítója.
  • while (/ N + / ig) {: ez meg fog felelni minden egymást követő N vagy n esetnek (a i zászló kis- és nagybetűket nem tesz szükségessé). A Perl az összes mérkőzés kezdő pozícióját a @ - speciális tömbben, a végállomásokat pedig a @ + -ban tárolja.
  • for (0 .. $ # -) {: az összes számot iterálja 0-tól a tömb végső indexéig ( $ # - ) @ - .
  • "$ F [0] \ t $ - [$ _] \ t $ + [$ _] \ n" nyomtatás: minimális ágyformátum-adatok nyomtatása. Az aktuális sorozat neve ( $ F [0] ), az aktuális mérkőzés kezdő pozíciója ( $ - [$ _] ) és a megfelelő véghelyzet ( $ + [$ _] ).

Éppen a fentieket futtattam a rendszeremen, és ez generálta ezt.

• @ user172818 kód nem számolhatja N-t a fa kezdetén
• @terdon-kód rossz ágyhelyzetet kap
• megváltoztatom a @terdon-kódot, mint az alábbiakban; és ez nálam működik

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