Come ottengo un file bed con l'elenco delle N-regioni di riferimento GRCh38? Queste sono le regioni in cui la sequenza è un tratto di Ns.
# se hai installato seqtk, salta le seguenti due righe git clone https://github.com/lh3/seqtkcd seqtk && rende # la riga di comando principale./ seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed L'opzione -n imposta la lunghezza minima di allungamento. Usa -p così com'è. È un po 'complicato spiegare cosa sta facendo.
Queste informazioni sono memorizzate nella rappresentazione del file a 2 bit della sequenza, quindi se ti capita di avere un file a 2 bit localmente (o vuoi scaricarne uno da UCSC) e hai installato py2bit (avrai bisogno della versione 3.0, dato che ho letteralmente appena aggiunto il supporto per questo):
import py2bittb = py2bit.open ("genome.2bit") of = open ("NNNNN.bed", "w") # Cambia i nomi dei file di input e output per chrom in tb.chroms (). keys (): blocks = tb.hardMaskedBlocks (chrom) per il blocco in blocchi: of.write ("{} \ t {} \ t {} \ n ".format (chrom, block [0], block [1])) of.close () tb.close () Se è utile, puoi usalo anche per ottenere tutte le regioni con maschera morbida. Sostituisci hardMaskedBlocks () con softMaskedBlocks () e assicurati di specificare storeMasked = True quando apri il file.
Usa twoBitInfo:
$ twoBitInfo file.fa -nBed output.bed Ad esempio, per ottenere tutte le regioni N-mascherate sul cromosoma Y (nota anche che puoi usare stdout come nome file per scrivere direttamente su stdout e usare l'URL come input, non è necessario scaricare il file a 2 bit):
$ twoBitInfo http: / /hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed stdout | head -3chrY 0 10000chrY 44821 94821chrY 133871 222346 Puoi scaricare twoBitInfo dalla directory appropriata al tuo sistema operativo qui: http://hgdownload.soe.ucsc.edu/admin / exe /
Potrei anche saltare. Ecco uno script Perl che ho scritto tempo fa per dividere una sequenza di fasta a tratti di Ns:
https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit. pl
L'ho appena modificato per sputare un file formattato BED che mostra le posizioni N su errore standard. Usalo come segue:
./fasta-nsplit.pl tig87_withN.fa 2>out.bed > out.split.fa Output (file BED) :
tig00000087 0 60tig00000087 900 960tig00000087 3840 3960tig00000087 14880 14940tig00000087 59520 59700tig00000087 93000 93120tig00000087 107880 107940tig00000087 109135 109140 : #! / usr / bin / perluse 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 + = lunghezza ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + length ($ nStretch)); $ cumLength + = lunghezza ($ nStretch); $ seq = $ newSeq; } printf (">% s \ n% s \ n", $ seqID, $ seq) if ($ seq); } $ seq = ""; $ shortSeqID = $ newShortID; $ seqID = $ newID; $ cumLength = 0; } altro {$ 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 + = lunghezza ($ seq); printf (STDERR "% s \ t% d \ t% d \ n", $ shortSeqID, $ cumLength, $ cumLength + length ($ nStretch)); $ cumLength + = lunghezza ($ nStretch); $ seq = $ newSeq; } printf (">% s \ n% s \ n", $ seqID, $ seq) if ($ seq);}
Ecco un modo per generarlo da soli dalla sequenza del genoma. Innanzitutto, converti il file fasta del genoma nel formato tbl ( <seq id> \ t<sequence>) , quindi usa perl per trovare le posizioni di inizio e fine di tutti i tratti consecutivi di N o n .
FastaToTbl hg38.fa | perl -F "\ t" -ane 'while (/ N + / ig) {for (0 .. $ # -) {print "$ F [0] \ t $ - [$ _] \ t $ + [$ _ ] \ n "}} '> hg38.n.bed FastaToTbl : questo è un script molto semplice che converte fasta in tbl. Salva semplicemente le righe sottostanti da qualche parte nel tuo $ PATH (ad esempio ~ / bin ) come FastaToTbl e rendilo eseguibile.
#! / 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 "} La magia del Perl.
-a fa sì che perl si comporti come awk , dividendo ogni riga di input sul valore fornito da -F (una scheda in questo caso) e salvare il risultato nel file aray @F . Quindi, $ F [0] sarà l'id della sequenza. while (/ N + / ig) {: corrisponderà a tutti i casi di N o n consecutivi (il i lo rende insensibile al maiuscolo / minuscolo). Perl memorizzerà le posizioni iniziali di tutte le corrispondenze nello speciale array @ - e le posizioni finali in @ + . for (0 .. $ # -) {: itera su tutti i numeri da 0 all'indice finale ( $ # - ) nell'array @ - . print "$ F [0] \ t $ - [$ _] \ t $ + [$ _] \ n" : stampa i dati sul formato minimo del letto. Il nome della sequenza corrente ( $ F [0] ), la posizione iniziale della corrispondenza corrente ( $ - [$ _] ) e la posizione finale corrispondente ( $ + [$ _] ). Ho appena eseguito quanto sopra sul mio sistema e ha generato questo.
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 "; }} "