Perl с извлечением последовательности FASTA имеет проблемы (только) с первой последовательностью

Question

Я использую функцию / подпрограмму extract_seq, доступную в Интернете, для извлечения последовательностей в файлах FASTA. Кратко:

• Последовательность начинается с первой строки, обозначенной '>', за которой следуют идентификатор и другая информация, разделенные пробелами
• Последующие строки (не начинающиеся с '>' имеют несколько строк
• Файл FASTA может иметь 1 или более последовательностей
• Ошибка в том, что в выводе есть дополнительный символ '>' для первой последовательности (только), что вызывает проблемы согласованности.

Программа отлично работает при извлечении последовательностей на основе идентификатора, за исключением дополнительного «>» в ​​случае первой последовательности. Не могли бы вы предложить решение, а также причину ошибки? Простое регулярное выражение решило бы проблему, но я не очень хорошо исправляю ошибки, которые не могу понять.

Скрипт Perl:

    #!/usr/bin/perl -w

    use strict;

    my $seq_all = "seq_all.fa";    # all proteins in fasta format

    foreach my $q_seq ("A0A1D8PC43","A0A1D8PC38") {
        print "Querying $q_seq\n";
        &extract_seq($seq_all, $q_seq);
    }

exit 0;

sub extract_seq
{
    open(my $fh, ">query.seq");

    my $seq_all = $_[0];
    my $lookup = $_[1];

    local $/ = "\n>";

    @ARGV = ($seq_all);
    while (my $seq = <>) {
        chomp $seq;
        my ($id) = $seq =~ /^>*(\S+)/;
        if ($id eq $lookup) {
            print "$seq\n";
            last;
        }
    }
}

Файл FASTA:

>A0A1D8PC43 A0A1D8PC43_CANAL Diphosphomevalonate decarboxylase
MYSASVTAPVNIATLKYWGKRDKSLNLPTNSSISVTLSQDDLRTLTTASASESFEKDQLW
LNGKLESLDTPRTQACLADLRKLRASIEQSPDTPKLSQMKLHIVSENNFPTAAGLASSAA
GFAALVSAIAKLYELPQDMSELSKIARKGSGSACRSLFGGFVAWEMGTLPDGQDSKAVEI
APLEHWPSLRAVILVVSDDKKDTPSTTGMQSTVATSDLFAHRIAEVVPQRFEAMKKAILD
KDFPKFAELTMKDSNSFHAVCLDSYPPIFYLNDTSKKIIKMVETINQQEVVAAYTFDAGP
NAVIYYDEANQDKVLSLLYKHFGHVPGWKTHYTAETPVAGVSRIIQTSIGPGPQETSESL
TK
>A0A1D8PC56 A0A1D8PC56_CANAL Uncharacterized protein OS=Candida
MSDTKKTTETDSEVGYLDIYLRFNDDMEKDYCFQVKTTTVFKDLYKVFRTLPISLRPSVF
YHAQPIGFKKSVSPGYLTQDGNFIFDEDSQKQAVPVNDNDLINETVWPGQLILPVWQFND
FGFYSFLAFLACWLYTDLPDFISPTPGICLTNQMTKLMAWVLVQFGKDRFAETLLADLYD
TVGVGAQCVFFGFHIIKCLFIFGFLYTGVFNPMRVFRLTPRSVKLDVTKEELVKLGWTGT
RKATIDEYKEYYREFKINQHGGMIQAHRAGLFNTLRNLGVQLESGEGYNTPLTEENKLRT
MRQIVEDAKKPDFKLKLSYEYFAELGYVFATNAENKEGSELAQLIKQYRRYGLLVSDQRI
KTVVRARKGETDEEKPKVEEVVEE
>A0A1D8PC67 A0A1D8PC67_CANAL Bfa1p OS=Candida albicans (strain
MVSDKLTLLRQFSEEDELFGDIEGIDYHDGETLKINKFSFPSSASSPSFAITGQSPNMRS
INGKRITRETLSEYSEENETDLTSEFSDQEFEWDGFNKNQSIYQQMNQRLIATKVAKQRE
AEREQRELMQKRHKDYDPNQTLRLKDFNKLTNENLTLLDQLDDEKTVNYEYVRDDVEDFA
QGFDKDFETKLRIQPSMPTLRSNAPTLKKYKSYGEFKCDNRVKQKLDRIPSFYNKNQLLS
KFKETKSYHPHHKKMGTVRCLNNNSEVPVTYPSISNMKLNKEKNRWEGNDIDLIRFEKPS
LITHKENKTKKRQGNMVYDEQNLRWINIESEHDVFDDIPDLAVKQLQSPVRGLSQFTQRT
TSTTATATAPSKNNETQHSDFEISRKLVDKFQKEQAKIEKKINHWFIDTTSEFNTDHYWE
IRKMIIEE
>A0A1D8PC38 A0A1D8PC38_CANAL Cta2p OS=Candida albicans (strain
MPENLQTRLHNSLDEILKSSGYIFEVIDQNRKQSNVITSPNNELIQKSITQSLNGEIQNF
HAILDQTVSKLNDAEWCLGVMVEKKKKHDELKVKEEAARKKREEEAKKKEEEAKKKAEEA
KKKEEEAKKAEEAKKAEEAKKVEEAAKKAEEAKKAEEEARKKAETAPQKFDNFDDFIGFD
INDNTNDEDMLSNMDYEDLKLDDKVPATTDNNLDMNNILENDESILDGLNMTLLDNGDHV
NEEFDVDSFLNQFGN

Edit: Проблема, как я уже объяснил выше, заключается в том, что в выводе есть дополнительный символ «>» для первой последовательности (только). Я не вижу причины для того же самого, и это вызывает много проблем. Выход:

Querying A0A1D8PC43
>A0A1D8PC43 A0A1D8PC43_CANAL Diphosphomevalonate decarboxylase
MYSASVTAPVNIATLKYWGKRDKSLNLPTNSSISVTLSQDDLRTLTTASASESFEKDQLW
LNGKLESLDTPRTQACLADLRKLRASIEQSPDTPKLSQMKLHIVSENNFPTAAGLASSAA
GFAALVSAIAKLYELPQDMSELSKIARKGSGSACRSLFGGFVAWEMGTLPDGQDSKAVEI
APLEHWPSLRAVILVVSDDKKDTPSTTGMQSTVATSDLFAHRIAEVVPQRFEAMKKAILD
KDFPKFAELTMKDSNSFHAVCLDSYPPIFYLNDTSKKIIKMVETINQQEVVAAYTFDAGP
NAVIYYDEANQDKVLSLLYKHFGHVPGWKTHYTAETPVAGVSRIIQTSIGPGPQETSESL
TK
Querying A0A1D8PC38
A0A1D8PC38 A0A1D8PC38_CANAL Cta2p OS=Candida albicans (strain
MPENLQTRLHNSLDEILKSSGYIFEVIDQNRKQSNVITSPNNELIQKSITQSLNGEIQNF
HAILDQTVSKLNDAEWCLGVMVEKKKKHDELKVKEEAARKKREEEAKKKEEEAKKKAEEA
KKKEEEAKKAEEAKKAEEAKKVEEAAKKAEEAKKAEEEARKKAETAPQKFDNFDDFIGFD
INDNTNDEDMLSNMDYEDLKLDDKVPATTDNNLDMNNILENDESILDGLNMTLLDNGDHV
NEEFDVDSFLNQFGN

Answer 1

Вы устанавливаете разделитель записей на \n>. Это не относится к первой последовательности.

Фиксированная кодовая последовательность:

...
chomp $seq;

# for first sequence
$seq =~ s/^>//;

my ($id) = $seq =~ /^(\S+)/;
if ($id eq $lookup) {
...

Обратите внимание, что ваша реализация крайне неэффективна, потому что она читает и анализирует содержимое файла для каждого запроса. Как насчет разделения загрузки / разбора и запросов на отдельные функции?

Альтернативное решение: дать полный список значений поиска загрузчику. Затем он будет заполнять массив ответов при обнаружении совпадений во время чтения файла.

Answer 2

$/ - разделитель входной записи, установка local $/="\n>"; приводит к тому, что вход разбивается на запись, заканчивающуюся \n>, после chomp окончание удаляется, однако />*(\S+)/ может не совпадать, поскольку > используется из предыдущей записи.

из FASTA wikipedia строка, начинающаяся с >, является комментарием и не всегда может быть идентификатором. Однако в случае, если это всегда так, можно исправить следующее.

my ($id,$seq) = $seq =~ /^>*(.*)\n(\S+)/;