Как извлечь идентификаторы из совпадающей строки в файлах FASTA, используя скрипт Perl?

Question

Я новичок в Perl программировании и застрял с моим сценарием.

Я пытаюсь найти мотив в файле FASTA и, если он найден, распечатать идентификатор содержащих его белков.

Я могу загрузить свой файл, но после размещения моего мотива ничего не происходит. Я получаю следующую ошибку: Use of uninitialized value $data[0] in concatenation (.) or string at test.pl line 36, <STDIN> line 2.

Вот мой код:

#!/usr/bin/perl -w
# Searching for motifs
print "Please type the filename of the protein sequence data: ";
$proteinfilename = <STDIN>;
# Remove the newline from the protein filename
chomp $proteinfilename;
# open the file, or exit
unless ( open(FA, $proteinfilename) ) {
print "Cannot open file \"$proteinfilename\"\n\n";
exit;
}
@protein = <FA>; # Read the protein sequence data from the file, and store it into the array variable @protein

my (@description, @ID, @data);

while (my $protein = <FA>) {
    chomp($protein);
    @description = split (/\s+/, $protein);
    push (@ID, $description[0]);   
}
# Close the file 
close FA;
my %params = map { $_ => 1 } @ID;
# Put the protein sequence data into a single string, as it's easier to search for a motif in a string than in an array of lines
$protein = join( '', @protein);
# Remove whitespace
$protein =~ s/\s//g;

# ask for a motif or exit if no motif is entered.
do {
print "Enter a motif to search for: ";
$motif = <STDIN>;
# Remove the newline at the end of $motif
chomp $motif;

# Look for the motif
@data = split (/\s+/, $protein);

if ( $protein =~ /$motif/ ) {
    print $description[0]."\n" if(exists($params{$data[0]}));
}
# exit on an empty user input
} until ( $motif =~ /^\s*$/ );
# exit the program
exit;

Пример ввода:

>sp|O60341|KDM1A_HUMAN Lysine-specific histone demethylase 1A OS=Homo sapiens OX=9606 GN=KDM1A PE=1 SV=2
MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP
RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY

Допустим, я хотел бы найти мотив PMET в заданной последовательности. Если он существует, я хочу получить ID в качестве вывода -> O60341

Answer 1

Если вы работаете с файлами fasta в perl, сделайте себе одолжение и используйте Bio Perl вместо того, чтобы пытаться свернуть свой код читателя. Делает вещи намного проще и чище.

#!/usr/bin/env perl
use strict;
use warnings; # Prefer this instead of the -w switch
use feature qw/say/;
use Bio::SeqIO;

my $reader = Bio::SeqIO->new(-format => "fasta",
                             # Use -file => "foo.fasta" to read from a file instead.
                             -fh => \*DATA);
my $motif = "PMET";

while (my $seq = $reader->next_seq) {
    if ($seq->seq =~ /\Q$motif/) {
        my $id = (split /\|/, $seq->primary_id)[1];
        say "Matched: $id";
    }
}

__DATA__
>sp|O60341|KDM1A_HUMAN Lysine-specific histone demethylase 1A OS=Homo sapiens OX=9606 GN=KDM1A PE=1 SV=2
MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP
RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY

Answer 2

После совпадения кода $motif и вывода $id

use strict;
use warnings;
use feature 'say';

my $motif = 'PMET';

while(<DATA>) {
    next unless /$motif/;
    my $id = (split '\|')[1];
    say "Found: $id";
}

__DATA__
sp|O60341|KDM1A_HUMAN Lysine-specific histone demethylase 1A OS=Homo sapiens OX=9606 GN=KDM1A PE=1 SV=2 MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY

Вывода

Found: O60341

Если вы хотите получить более подробную информацию, следующий код выполняет эту функцию

use strict;
use warnings;
use feature 'say';

use Data::Dumper;

my $debug = 1;

my $motif = 'PMET';

while(<DATA>) {
    next unless /$motif/;
    my $record = fill_record($_);
    say 'Found: ' . $record->{id};
    say Dumper($record) if $debug;
}

sub fill_record {
    my $data = shift;

    my(%record,@dna);

    (@record{qw/x id/}, $data) = split '\|', $data;
    ($data,@dna)  = $data =~ /(.*?) (\S+) (\S+)$/;
    $record{dna}  = \@dna;
    $data =~ /(.*?) \w+=/;
    $record{desc} = $1;

    %{$record{group}} = $data =~ /\b(\w+)=([^=]+) (?=\w+=)/g;

    return \%record;
}

__DATA__
sp|O60341|KDM1A_HUMAN Lysine-specific histone demethylase 1A OS=Homo sapiens OX=9606 GN=KDM1A PE=1 SV=2 MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY

Выход

Found: O60341
$VAR1 = {
          'x' => 'sp',
          'id' => 'O60341',
          'group' => {
                       'OX' => '9606',
                       'OS' => 'Homo sapiens',
                       'GN' => 'KDM1A',
                       'PE' => '1'
                     },
          'desc' => 'KDM1A_HUMAN Lysine-specific histone demethylase 1A',
          'dna' => [
                     'MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP',
                     'RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY'
                   ]
        };

Answer 3

Я кодировал здесь для однострочного входного файла для примеров.

my $motif = <STDIN>;
chomp($motif);

my $str = "sp|O60341|KDM1A_HUMAN Lysine-specific histone demethylase 1A OS=Homo sapiens OX=9606 GN=KDM1A PE=1 SV=2 MLSGKKAAAAAAAAAAAATGTEAGPGTAGGSENGSEVAAQPAGLSGPAEVGPGAVGERTP RKKEPPRASPPGGLAEPPGSAGPQAGPTVVPGSATPMETGIAETPEGRRTSRRKRAKVEY";

if($str=~m/$motif/)
{
    if($str=~m/^([^|]+)\|([^|]+)\|/gm)
    {
        print "Expected Value: $2\n";
    }
}
else { print "Not matched...\n";  }

Input>$: PMET

Expected Value: O60341