Длина последовательности файла FASTA

У меня есть следующий файл FASTA:

>header1
CGCTCTCTCCATCTCTCTACCCTCTCCCTCTCTCTCGGATAGCTAGCTCTTCTTCCTCCT
TCCTCCGTTTGGATCAGACGAGAGGGTATGTAGTGGTGCACCACGAGTTGGTGAAGC
>header2
GGT
>header3
TTATGAT

Мой желаемый результат:

>header1
117
>header2
3
>header3
7
# 3 sequences, total length 127.

Это мой код:

awk '/^>/ {print; next; } { seqlen = length($0); print seqlen}' file.fa

Вывод, который я получаю с этим кодом:

>header1
60
57
>header2
3
>header3
7

Мне нужна небольшая модификация, чтобы иметь дело с несколькими строками последовательности.

Мне также нужен способ получить общие последовательности и общую длину. Любое предложение будет приветствоваться... В bash или awk, пожалуйста. Я знаю, что это легко сделать в Perl/BioPerl, и на самом деле у меня есть скрипт, чтобы сделать это таким образом.


bash awk fasta
person cucurbit    schedule 02.06.2014    source источник

Ответы (3)


arrow_upward
19
arrow_downward

Решение awk / gawk может состоять из трех этапов:

  1. Каждый раз при обнаружении header следует выполнять следующие действия:

    • Print previous seqlen if exists.
    • Печать тега.
    • Инициализировать последовательность.
  2. Для sequence строк нам просто нужно накопить итоги.
  3. Наконец, на этапе END мы печатаем остаточную последовательность.

Код с комментариями:

awk '/^>/ { # header pattern detected
        if (seqlen){
         # print previous seqlen if exists 
         print seqlen
         }

         # pring the tag 
         print

         # initialize sequence
         seqlen = 0

         # skip further processing
         next
      }

# accumulate sequence length
{
seqlen += length($0)
}
# remnant seqlen if exists
END{if(seqlen){print seqlen}}' file.fa

Однострочник:

awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' file.fa

Для итогов:

awk '/^>/ { if (seqlen) {
              print seqlen
              }
            print

            seqtotal+=seqlen
            seqlen=0
            seq+=1
            next
            }
    {
    seqlen += length($0)
    }     
    END{print seqlen
        print seq" sequences, total length " seqtotal+seqlen
    }' file.fa
person Juan Diego Godoy Robles    schedule 02.06.2014
comment
Да, это работает. Но также мне нужна последняя строка с общим количеством последовательностей и общей длиной, следуя примеру: 3 последовательности, общая длина 127. (Извините, в вопросе это за # ) - person cucurbit; 02.06.2014
comment
Просто незначительное изменение форматирования. awk '/^>/ {if (seqlen) print seqlen;print;seqlen=0;next} {seqlen+=length($0)}END{print seqlen}' - person Jotne; 02.06.2014
comment
Очень старая тема, но это как раз то, что мне нужно! Можно ли также настроить команду awk с одним вкладышем, чтобы она печатала количество нуклеотидов в той же строке, что и заголовок, а не в новой строке, например: ›header1 60 - person Gravel; 17.05.2017
comment
Конечно, @Gravel просто поместите ; printf $0" " ; вместо ; print ; - person Juan Diego Godoy Robles; 22.05.2017
comment
это будет включать символы новой строки в счет - person Brian Wiley; 08.09.2020
comment
извините, может быть, это подсчет возврата каретки. просто вставьте BEGIN{RS="\r\n"} в начале - person Brian Wiley; 08.09.2020

arrow_upward
1
arrow_downward

Быстрый способ с любым awk:

awk '/^>/{if (l!="") print l; print; l=0; next}{l+=length($0)}END{print l}' file.fasta

Вас также может заинтересовать BioAwk, это адаптированная версия awk, настроенная для обработки файлов FASTA.

bioawk -c fastx '{print ">" $name ORS length($seq)}' file.fasta

Примечание. BioAwk основан на awk Брайана Кернигана, задокументированный в "Язык программирования AWK", Аль Ахо, Брайан Керниган и Питер Вайнбергер (Addison-Wesley, 1988, ISBN 0-201-07981 -Х). Я не уверен, совместима ли эта версия с POSIX. .

person kvantour    schedule 12.12.2018

arrow_upward
0
arrow_downward

Я хотел поделиться некоторыми изменениями в ответе klashxx, которые могут быть полезны. Его вывод отличается тем, что он печатает идентификатор последовательности и ее длину в одной строке. Он больше не является однострочным, поэтому недостатком является то, что вам придется сохранять его как файл сценария.

Он также анализирует идентификатор последовательности из строки заголовка на основе пробелов (chrM в >chrM gi|251831106|ref|NC_012920.1|). Затем вы можете выбрать конкретную последовательность на основе идентификатора, установив переменную target следующим образом: $ awk -f seqlen.awk -v target=chrM seq.fa.

BEGIN {
  OFS = "\t"; # tab-delimited output
}
# Use substr instead of regex to match a starting ">"
substr($0, 1, 1) == ">" {
  if (seqlen) {
    # Only print info for this sequence if no target was given
    # or its id matches the target.
    if (! target || id == target) {
      print id, seqlen;
    }
  }
  # Get sequence id:
  # 1. Split header on whitespace (fields[1] is now ">id")
  split($0, fields);
  # 2. Get portion of first field after the starting ">"
  id = substr(fields[1], 2);
  seqlen = 0;
  next;
}
{
  seqlen = seqlen + length($0);
}
END {
  if (! target || id == target) {
    print id, seqlen;
  }
}
person Nick S    schedule 16.02.2015