Длина последовательности файла 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, и на самом деле у меня есть скрипт, который делает это таким образом.

Ответ 1

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

При каждом обнаружении header эти действия должны выполняться:
- Печать предыдущего seqlen , если существует.
- Печатать тег.
- Инициализировать seqlen.
Для строк sequence нам просто нужно аккумулировать итоговые значения.
Наконец, на этапе END мы печатаем остатки seqlen.

Прокомментированный код:

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

A oneliner:

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

Ответ 2

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

Он также анализирует идентификатор последовательности из строки заголовка на основе пробелов (chrM in >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;
  }
}

Ответ 3

Быстрый способ с любым 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-X). Я не уверен, что эта версия совместима с POSIX.