Liczenie określonego kolejnego znaku z jego pozycją i długością występowania

3

Mam plik sekwencji i chcę policzyć kolejny znak „N” z jego pozycją wystąpienia i długością Powiedzmy, że mam plik o nazwie mySequence.fastatakiej:

>sequence-1
ATCGCTAGCATNNNNNNNNNNNNNNCTAGCATCATGCNNNNNNATACGCATCACANNNNNNNNNCgcatATCAC

i oczekiwany wynik powinien wyglądać następująco:

Position 12 N 14
Position 38 N 6
Position 56 N 9

Uprzejmie pomóż mi rozwiązać ten problem awklub sedpodając nazwę mojego plikumySequence.fasta

awk sed perl bioinformatics Początkujący bioinformatyk
źródło
1
Czy plik Fasta zawiera tylko sekwencje o długości jednej linii, czy też rozwiązanie musi obejmować przypadek, w którym ciąg Ns owija się wokół końca linii? Ponadto, czy w pliku jest jedna lub kilka sekwencji?
Kusalananda
Moim celem jest tylko jedna sekwencja w pliku. Ponadto prawdą jest, że istnieje sekwencja wielu fasta (jeden plik zawierający wiele sekwencji zaczynających się na literę „>”. Więc kolejnym możliwym rozwiązaniem dla tego samego zadania byłoby dodanie identyfikatora sekwencji na wyjściu.
Budding-bioinformatician

Odpowiedzi:

10

Możesz to zrobić za pomocą awk, którego match()ustawienie RSTARTi RLENGTHzmienna jest do tego całkiem użyteczne:

<mySequence.fasta awk -v C=N '{
  i=0
  while (match($0, C "+")) {
    printf "Position %d %s %d\n", i+RSTART, C, RLENGTH
    i += RSTART+RLENGTH-1
    $0 = substr($0, RSTART+RLENGTH)
  }}'

Lub przy perlużyciu tablic @-i, @+które rejestrują początek i koniec meczów:

perl -ne 'printf "Position %d N %d\n", $-[0]+1, $+[0]-$-[0] while /N+/g'

Kolejne nieco szybsze (przynajmniej w mojej wersji perl) perlpodejście z wykorzystaniem ( eksperymentalnego ) (?{...})operatora wyrażenia regularnego:

perl -ne '0 while /N(?{$s=pos})N*(?{printf "Position %d N %s\n", $s, pos()-$s+1})/g'
Stéphane Chazelas
źródło
Rozwiązanie Perla jest szybsze niż awk
Budding-bioinformatician
2
@ Budding-bioinformatician, dla gawk, zmień ustawienia regionalne na C ( LC_ALL=C awk...), aby przyspieszyć (po czym uważam, że wydajność jest porównywalna perl). Uważam, że mawkjest to znacznie szybsze niż gawklub perlna tym. Możesz przyspieszyć to jeszcze bardziej, jeśli na stałe kodujesz Nw kodzie jak w perl.
Stéphane Chazelas
2

Inne awkrozwiązanie:

awk -F '' '{for(i=1;i<=NF;i++){ if($i=="N"&&!sPOS) sPOS=i;
   if (i==NF &&sPOS && $NF=="N"){LN++}; if($i=="N" &&sPOS && i<NF) {LN++}
   else if(sPOS) {printf("Position %d N %d\n", sPOS, LN); LN=sPOS=0} }
}' infile.txt

Ponieważ cała awkimplementacja nie obsługuje pustego FS ( -F ''), poniżej znajduje się poprawiony skrypt, który jest zgodny:

awk -F'N' '{sPOS=0;for(i=2;i<=NF;i++){ if($i==""&&!sPOS) sPOS=(i-1)+length($(i-1));
    if($i=="" &&sPOS && NF!=i) {LN++} 
    else if(sPOS) {printf("Position %d N %d\n", sPOS, ++LN); sPOS+=LN+length($i); LN=0} }
}' infile.txt

Przykładowe dane wejściowe:

>sequence-1
ATCGCTAGCATNNNNNNNNNNNNNNCTAGCATCATGCNNNNNNATACGCATCACANNNNNNNNNCgcatATCACNN
N
AN
NNA

Wynik to:

Position 12 N 14
Position 38 N 6
Position 56 N 9
Position 75 N 2
Position 1 N 1
Position 2 N 1
Position 1 N 2
αғsнιη
źródło