Обрежьте первые N баз в мультифайле fasta с помощью awk и распечатайте в формате максимальной ширины

Фон

Формат multi fasta содержит несколько записей последовательностей, каждая запись начинается с однострочного описания, за которым следуют несколько строк последовательности (РНК, ДНК, белок). Строка описания имеет в начале символ больше, после которого следует ">" - идентификатор последовательности, а оставшаяся часть строки содержит описание записи (оба значения необязательны).

в файлах fasta обычно последовательности строк должны быть отформатированы до максимальной ширины

пример ввода с максимальной шириной = "70":

>gi|304322925|ref|YP_003856771.1| NADH [Lynx rufus]
MMTYIVFILSTIFVVSFVGFSSKPSPIYGGFGLIVAGGIGCGIVLNFGGSFLGLMVFLIYLGGMLVVFGY
TTAMATEPYPEAWTSNKAVLGMLITGILAELLTACYILKEDEIEVVFKFNGAGDWVIYDTGDSGFFSEEA
MGIAALYSYGTWLVVVTGWSLLIGVLVIMEVTRGN
>gi|295065592|ref|YP_003587393.1| NADH [Nomascus siki]
MTYTLFLLSVILVMGFVGFSSKPSPIYGGLVLVVSGVVGCAVILNCGGGYLGLMVFLIYLGGMMVVFGYT
TAMAIEEYPEAWGSGVEVLVGVLVGLAMEVGLVLWAKECDGLVMVLNFNNMGSWVIYEGEGSGLIREDSI
GAGALYDYGRWLVVVTGWTLLVGVYIVIEIARGN
>gi|295065550|ref|YP_003587316.1| NADH (mitochondrion) [Symphalangus syndactylus]
MTYTLFLLSVILVMGFVGFSSKPSPIYGGLVLVVSGVVGCAIILDCGGGYLGLMVFLIYLGGMMVVFGYT
TAMAIEEYPEAWGSGVEVLVGVLVGLAMEVGLVLWAKEYDGLVVVLNFNNMGSWVIYEGEGSGLIREDSI
GAGALYDYGRWLVVVTGWTLLVGVYIVIEIARGN

Мои усилия

Я написал gawk-скрипт, он оставляет обрезку в последовательности (обрезает первые N оснований), и если длина последовательности меньше, чем обрезка, то не печатает ее.

gawk -v left_trim=15 -v width=70 '
  BEGIN{RS=">"}
  NR==1{next};        #for blank record
  {
    split ($0, a, "\n"); 
    sequence=""; 
    for(i=2; i<=length(a); i++){
      sequence=sequence a[i];
    }; 
    if(length(sequence)>left_trim) {
      printf(">%s\n%s\n",a[1], substr(sequence, left_trim+1));
    }
  }' test.fasta

результат

>gi|304322925|ref|YP_003856771.1| NADH [Lynx rufus]
SFVGFSSKPSPIYGGFGLIVAGGIGCGIVLNFGGSFLGLMVFLIYLGGMLVVFGYTTAMATEPYPEAWTSNKAVLGMLITGILAELLTACYILKEDEIEVVFKFNGAGDWVIYDTGDSGFFSEEAMGIAALYSYGTWLVVVTGWSLLIGVLVIMEVTRGN
>gi|295065592|ref|YP_003587393.1| NADH [Nomascus siki]
FVGFSSKPSPIYGGLVLVVSGVVGCAVILNCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSGVEVLVGVLVGLAMEVGLVLWAKECDGLVMVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVVTGWTLLVGVYIVIEIARGN
>gi|295065550|ref|YP_003587316.1| NADH (mitochondrion) [Symphalangus syndactylus]
FVGFSSKPSPIYGGLVLVVSGVVGCAIILDCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSGVEVLVGVLVGLAMEVGLVLWAKEYDGLVVVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVVTGWTLLVGVYIVIEIARGN

ожидается

>gi|304322925|ref|YP_003856771.1| NADH [Lynx rufus]
SFVGFSSKPSPIYGGFGLIVAGGIGCGIVLNFGGSFLGLMVFLIYLGGMLVVFGYTTAMATEPYPEAWTS
NKAVLGMLITGILAELLTACYILKEDEIEVVFKFNGAGDWVIYDTGDSGFFSEEAMGIAALYSYGTWLVV
VTGWSLLIGVLVIMEVTRGN
>gi|295065592|ref|YP_003587393.1| NADH [Nomascus siki]
FVGFSSKPSPIYGGLVLVVSGVVGCAVILNCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKECDGLVMVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN
>gi|295065550|ref|YP_003587316.1| NADH (mitochondrion) [Symphalangus syndactylus]
FVGFSSKPSPIYGGLVLVVSGVVGCAIILDCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKEYDGLVVVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN

Вопросы

  • как я могу реализовать формат ширины? я пытаюсь с ">%s\n%.70s\n", напечатать первые 70 букв, ">%s\n%70s\n" напечатать все

  • как я могу улучшить это соединение? выходит из функции?

    sequence=""; 
for(i=2; i<=length(a); i++){
  sequence=sequence a[i];
};

awk fasta gawk
person Jose Ricardo Bustos M.    schedule 15.05.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Чтобы ответить на ваши конкретные вопросы, вы можете указать ширину поля вывода, используя модификатор формата *:

$ awk 'BEGIN{printf "%s\n", "foo"}'
foo
$ awk 'BEGIN{printf "%*s\n", 10, "foo"}'
       foo

и нет, нет функции join для объединения массивов в строку (противоположность split()), НО если вы просто позволите awk разделить запись на поля вместо того, чтобы вручную разбивать запись на массив элементов, тогда вы можете просто присвойте новое значение любому полю, и awk перекомпилирует поля в $0, как я делаю в качестве преднамеренного побочного эффекта в моем первом решении ниже, когда я удаляю первое поле/строку с $1 = "".

Вот как я бы подошел к задаче в любом случае:

$ cat tst.awk
BEGIN { RS=">"; FS="\n"; OFS="" }
NR>1 {
    print RS $1
    $1 = ""
    for ( start=left_trim+1; start<=length(); start+=width ) {
        print substr($0,start,width)
    }
}

$ awk -v left_trim=15 -v width=70 -f tst.awk file
>gi|304322925|ref|YP_003856771.1| NADH [Lynx rufus]
SFVGFSSKPSPIYGGFGLIVAGGIGCGIVLNFGGSFLGLMVFLIYLGGMLVVFGYTTAMATEPYPEAWTS
NKAVLGMLITGILAELLTACYILKEDEIEVVFKFNGAGDWVIYDTGDSGFFSEEAMGIAALYSYGTWLVV
VTGWSLLIGVLVIMEVTRGN
>gi|295065592|ref|YP_003587393.1| NADH [Nomascus siki]
FVGFSSKPSPIYGGLVLVVSGVVGCAVILNCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKECDGLVMVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN
>gi|295065550|ref|YP_003587316.1| NADH (mitochondrion) [Symphalangus syndactylus]
FVGFSSKPSPIYGGLVLVVSGVVGCAIILDCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKEYDGLVVVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN

или, если вы предпочитаете:

$ cat tst.awk
/^>/ { prtRec(); rec=""; print; next }
{ rec = rec $0 }
END { prtRec() }
function prtRec(        start) {
    for ( start=left_trim+1; start<=length(rec); start+=width ) {
        print substr(rec,start,width)
    }
}

$ awk -v left_trim=15 -v width=70 -f tst.awk file
>gi|304322925|ref|YP_003856771.1| NADH [Lynx rufus]
SFVGFSSKPSPIYGGFGLIVAGGIGCGIVLNFGGSFLGLMVFLIYLGGMLVVFGYTTAMATEPYPEAWTS
NKAVLGMLITGILAELLTACYILKEDEIEVVFKFNGAGDWVIYDTGDSGFFSEEAMGIAALYSYGTWLVV
VTGWSLLIGVLVIMEVTRGN
>gi|295065592|ref|YP_003587393.1| NADH [Nomascus siki]
FVGFSSKPSPIYGGLVLVVSGVVGCAVILNCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKECDGLVMVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN
>gi|295065550|ref|YP_003587316.1| NADH (mitochondrion) [Symphalangus syndactylus]
FVGFSSKPSPIYGGLVLVVSGVVGCAIILDCGGGYLGLMVFLIYLGGMMVVFGYTTAMAIEEYPEAWGSG
VEVLVGVLVGLAMEVGLVLWAKEYDGLVVVLNFNNMGSWVIYEGEGSGLIREDSIGAGALYDYGRWLVVV
TGWTLLVGVYIVIEIARGN
person Ed Morton    schedule 15.05.2015
comment
$1="" лучше sub(/[^\n]+\n/,"") только из-за OFS="". Без этого присваивание привело бы к появлению нежелательного пустого символа в начале строки. Комбинация $1="" и OFS="" делает gsub() ненужным - $1="" говорит awk перекомпилировать $0, а OFS="" говорит ему просто удалить символы новой строки (из FS="\n") между полями во время этой перекомпиляции вместо замены их каким-либо другим символом (пустым по умолчанию). - person Ed Morton; 15.05.2015