从包含其他信息的巨大文本文件中定位和获取序列答案

【问题标题】：Locating and obtaining a sequence from a huge text file which includes other information从包含其他信息的巨大文本文件中定位和获取序列
【发布时间】：2014-03-15 09:11:19
【问题描述】：

我得到了 2 个 P3DB 格式的文本文件（第一个称为“位点”具有蛋白质 ID，第二个称为“蛋白质”具有相应蛋白质 ID 的序列）。我必须使用 Perl 脚本将这两个转换为单个 PhosphoSitePlus 文本文件。

我已经弄清楚如何转换大部分信息。现在我必须在“蛋白质”文件中找到序列和生物体名称。我有来自“站点”文件的某个蛋白质 ID 号（例如，2329）。现在我必须在“蛋白质”文件中搜索这个数字。我找到了这个数字，但在它下面有很多不必要的数据，在某个地方是有机体名称后面的序列。我不确定如何获取序列，因为我不知道它什么时候开始或停止。

有一种模式，所有序列都以制表符空格和“M”开头。但是，最后一个残基可以是任何东西。此外，有机体名称可能紧跟在序列之后（不带空格）或在序列之后的多个空格之后。

我希望能够保存与蛋白质 ID 对应的完整序列，以便能够在该序列中的给定位置（编号）找到单个残基。或者能够找到一个残基而无需先保存序列。

这大概是我想出的。

#!/usr/bin/perl
use strict;
use warnings;

my $filename = $ARGV[0];
my $filename2 = $ARGV[1];

open AFILE, "$filename"; #open file - sites file 
open BFILE, "$filename2"; #open second file - proteins file 
open NEWFILE, ">PhosphoSitePlus.txt"; #make a new file to save

my @chunks = ();
my @blines = <BFILE>;
my $sequence;
my $res;
my $organism;
my $PID;
my $ACC;
my $Psite;

print NEWFILE "Accession        Modified Residue        Site Group ID   Organism        Sequence";

while (defined(my $line = <AFILE>)) { #iterate through lines of the _sites file
    my @chunks = split ' ', $line; #split columns into an array
    next if ( $chunks[0] =~  /P3DB/ ); #skip the first line (the headers)
    $PID = $chunks[0]; #save the Protein ID
    $ACC = $chunks[1]; #save the Accession number
    $Psite = $chunks[3]; #save the site that is phosphorylated

    foreach (0 .. $#blines){ #iterate through p3db_proteins

        my @b = split ' ', $blines[$_]; #split columns of lines into array
        if ($b[0] =~/^$PID$/){  #if the first column is equal to the protein ID the sequence is under
            next if ($blines[$_] !~ /^\s+M\w{20,}$/);
            if ( $blines[$_] =~ /^\tM\w{20,}$/ ){ #the start of the sequence(tab and an M, followed by 20 or more characters?)
                #save this sequence, or find the residue that is at the position of $Psite
            }
        }
    }
}

以下是蛋白质文件中序列的一些示例（2329 和 2330 是蛋白质 ID，以 M 开头的是序列和生物名称）。

2329    EMBL:AAM13013.1;EMBL:AAM65937.1;EMBL:AAP13391.1;EMBL:AEC10514.1;Ensembl Genomes:AT2G45140;TAIR:At2g45140;TAIR:AT2G45140.1;Ensembl Genomes:AT2G45140.1;TAIR:AT2G45140.1;PIR:H84886;IPI:IPI00531520.1;Refseq:NM_130077.2;Refseq:NP_182039.1;Swissport:Q9SHC8.1;UniParc:UPI00000A0803;Swissport:VAP12_ARATH    plant VAP homolog 12    MSNELLTIDPVDLQFPFELKKQISCSLYLGNKTDNYVAFKVKTTNPKKYCVRPNTGVVHPRSSSEVLVTMQAQKEAPADLQCKDKFLLQCVVASPGATPKDVTHEMFSKEAGHRVEETKLRVVYVAPPRPPSPVREGSEEGSSPRASVSDNGNASDFTAAPRFSADRVDAQDNSSEARALVTKLTEEKNSAVQLNNRLQQELDQLRRESKRSKSGGIPFMYVLLVGLIGLILGYIMKRT Arabidopsis thaliana    19376835;19253305;19245862;18463617;17651370;17317660;15308754
2330    EMBL:AEE76598.1;Ensembl Genomes:AT3G22180;TAIR:At3g22180;TAIR:AT3G22180.1;Ensembl Genomes:AT3G22180.1;TAIR:AT3G22180.1;EMBL:BAB03066.1;IPI:IPI00547221.2;Refseq:NM_113115.3;Refseq:NP_188857.1;Swissport:Q9LIE4.2;UniParc:UPI00001634CF;Swissport:ZDHC8_ARATH   DHHC-type zinc finger family protein    MVRKHGWQLPAHTLQVIAITVFCLLVVAFYAFFAPFVGGRIWEYVLIGVYSPVAILVFVLYVRCTAINPADPRIMSIFDTGVNGDGMVRGLSRNYDETGSQLQASPSVVSRSSTVAGNSSVKGSVEDAQRVESVSRRSCYNPLAVFCYVFVVEDCRKKEGPAEEQGNSEEALFCTLCNCEVRKFSKHCRSCDKCVDCFDHHCKWLNNCVGRKNYVTFVSLMSASLLWLIIEAAVGIAVIVRVFVNKQTMETEIVNRLGNSFSRAPLAAVVGLCTAVAIFACFPLGELLFFHMLLIKKGITTYEYVVAMRAMSEAPDGASVDEEIQNVLYSPTGSATTGFSGGSSLGLPYRGVWCTPPRVFDNQDEVIPHLDPCMVPSTVDPDAPGSEKGTKALKRPVKRNAWKLAKLDPNEAARAAARARASSSVLRPIDNRHLPDNDLSSIGTVSIISSVSTDANVAASKEIRNNDLRSSLSRNSFAPSQGSRDEYDTGSHGMSNLSSPSHVHESVTLAPLPQNPTIVGNRFTATSHHMHSTFDDKVLHRGNDADPLFLFAPATSHLRDVRKTSVVWDPEAGRYVSAPVTTTSEVRNRLLNPSSQTASTQNPRPILPAHDSSSGSSALRDPLPLHQAERRLTYTGDSIFYGGPLINIPTRDTPRSGRGLVRDVQDRLASTVHRDARIRRDSTSNQLPVFAPGGLGANSQTGSNIK  Arabidopsis thaliana    17317660

【问题讨论】：

请您编辑您的问题以便正确显示数据？您使用了直角括号>，用于引用散文部分。它将多个空格减少为单个空格，并将行一起运行而不会中断以及使用可变宽度字体，因此很难看到原始数据的样子。您可以通过将块的每一行缩进四个空格来标记一段代码或数据。如果您直接从文件中粘贴文本，选择它，然后按CTRL-K，可能会更简单。这将为您缩进整个块。
我已经使用 CTRL-K 编辑了数据。谢谢。
@Borodin 感谢CTRL-K 的把戏。

标签： perl

【解决方案1】：

您正在使用匹配行首的^ 和匹配行尾的$，这意味着该行上唯一的内容是您要查找的以 M 开头的序列。

但如果你的数据如图所示，那就不是这样了。

另外，如果你的文件和你说的一样大，我建议不要把整个文件都塞进内存，而是一行一行地做。

下面的正则表达式应该做你想做的事。

#!/usr/bin/perl

use strict;
use warnings;
use Data::Dumper;

my $id = 2329;
my @info;
while(my $line = <DATA>){
    if($line =~ /(2329|2330)/){
        push @info, $line =~ /\s+(M\w+)/;
    }
}

print Dumper @info;

__DATA__
2329    EMBL:AAM13013.1;EMBL:AAM65937.1;EMBL:AAP13391.1;EMBL:AEC10514.1;Ensembl Genomes:AT2G45140;TAIR:At2g45140;TAIR:AT2G45140.1;Ensembl Genomes:AT2G45140.1;TAIR:AT2G45140.1;PIR:H84886;IPI:IPI00531520.1;Refseq:NM_130077.2;Refseq:NP_182039.1;Swissport:Q9SHC8.1;UniParc:UPI00000A0803;Swissport:VAP12_ARATH    plant VAP homolog 12    MSNELLTIDPVDLQFPFELKKQISCSLYLGNKTDNYVAFKVKTTNPKKYCVRPNTGVVHPRSSSEVLVTMQAQKEAPADLQCKDKFLLQCVVASPGATPKDVTHEMFSKEAGHRVEETKLRVVYVAPPRPPSPVREGSEEGSSPRASVSDNGNASDFTAAPRFSADRVDAQDNSSEARALVTKLTEEKNSAVQLNNRLQQELDQLRRESKRSKSGGIPFMYVLLVGLIGLILGYIMKRT Arabidopsis thaliana    19376835;19253305;19245862;18463617;17651370;17317660;15308754
2330    EMBL:AEE76598.1;Ensembl Genomes:AT3G22180;TAIR:At3g22180;TAIR:AT3G22180.1;Ensembl Genomes:AT3G22180.1;TAIR:AT3G22180.1;EMBL:BAB03066.1;IPI:IPI00547221.2;Refseq:NM_113115.3;Refseq:NP_188857.1;Swissport:Q9LIE4.2;UniParc:UPI00001634CF;Swissport:ZDHC8_ARATH   DHHC-type zinc finger family protein    MVRKHGWQLPAHTLQVIAITVFCLLVVAFYAFFAPFVGGRIWEYVLIGVYSPVAILVFVLYVRCTAINPADPRIMSIFDTGVNGDGMVRGLSRNYDETGSQLQASPSVVSRSSTVAGNSSVKGSVEDAQRVESVSRRSCYNPLAVFCYVFVVEDCRKKEGPAEEQGNSEEALFCTLCNCEVRKFSKHCRSCDKCVDCFDHHCKWLNNCVGRKNYVTFVSLMSASLLWLIIEAAVGIAVIVRVFVNKQTMETEIVNRLGNSFSRAPLAAVVGLCTAVAIFACFPLGELLFFHMLLIKKGITTYEYVVAMRAMSEAPDGASVDEEIQNVLYSPTGSATTGFSGGSSLGLPYRGVWCTPPRVFDNQDEVIPHLDPCMVPSTVDPDAPGSEKGTKALKRPVKRNAWKLAKLDPNEAARAAARARASSSVLRPIDNRHLPDNDLSSIGTVSIISSVSTDANVAASKEIRNNDLRSSLSRNSFAPSQGSRDEYDTGSHGMSNLSSPSHVHESVTLAPLPQNPTIVGNRFTATSHHMHSTFDDKVLHRGNDADPLFLFAPATSHLRDVRKTSVVWDPEAGRYVSAPVTTTSEVRNRLLNPSSQTASTQNPRPILPAHDSSSGSSALRDPLPLHQAERRLTYTGDSIFYGGPLINIPTRDTPRSGRGLVRDVQDRLASTVHRDARIRRDSTSNQLPVFAPGGLGANSQTGSNIK  Arabidopsis thaliana    17317660

另外，在您的while(defined(my $line = <AFILE>)) 行上，perl 已经知道您在这里的意思，您不需要包含已定义的部分。如果它没有收到数据/或到达文件末尾，$line 将评估为 false。

【讨论】：

【解决方案2】：

好的，我想我明白你想要什么了。

您显示的少量数据似乎表现得相当好，因为两行都有六个制表符分隔的字段。 ID 是第一个字段，序列是第四个。您还提到了有机体名称，我认为这是第五个字段，但您没有明确说明您是否对这个值感兴趣。

既然您说您的数据集很大，要从中获得任何类型的性能，您需要从蛋白质文件中构建一个 hash，其中蛋白质 ID 作为键，序列和有机体名称作为值。

此后，您可以只查找给定蛋白质 ID 的数据，而不是搜索文件。

这个简短的程序展示了一个例子

use strict;
use warnings;

my %proteins;

open my $proteins_fh, '<', 'proteins.pdb' or die $!;

while (<$proteins_fh>) {
  my @fields = split /\t/;
  my ($pid, $sequence) = @fields[0, 3];
  $proteins{$pid} = $sequence;
}

use Data::Dump;
dd \%proteins;

输出

{
  2329 => [
            "MSNELLTIDPVDLQFPFELKKQISCSLYLGNKTDNYVAFKVKTTNPKKYCVRPNTGVVHPRSSSEVLVTMQAQKEAPADLQCKDKFLLQCVVASPGATPKDVTHEMFSKEAGHRVEETKLRVVYVAPPRPPSPVREGSEEGSSPRASVSDNGNASDFTAAPRFSADRVDAQDNSSEARALVTKLTEEKNSAVQLNNRLQQELDQLRRESKRSKSGGIPFMYVLLVGLIGLILGYIMKRT",
            "Arabidopsis thaliana",
          ],
  2330 => [
            "MVRKHGWQLPAHTLQVIAITVFCLLVVAFYAFFAPFVGGRIWEYVLIGVYSPVAILVFVLYVRCTAINPADPRIMSIFDTGVNGDGMVRGLSRNYDETGSQLQASPSVVSRSSTVAGNSSVKGSVEDAQRVESVSRRSCYNPLAVFCYVFVVEDCRKKEGPAEEQGNSEEALFCTLCNCEVRKFSKHCRSCDKCVDCFDHHCKWLNNCVGRKNYVTFVSLMSASLLWLIIEAAVGIAVIVRVFVNKQTMETEIVNRLGNSFSRAPLAAVVGLCTAVAIFACFPLGELLFFHMLLIKKGITTYEYVVAMRAMSEAPDGASVDEEIQNVLYSPTGSATTGFSGGSSLGLPYRGVWCTPPRVFDNQDEVIPHLDPCMVPSTVDPDAPGSEKGTKALKRPVKRNAWKLAKLDPNEAARAAARARASSSVLRPIDNRHLPDNDLSSIGTVSIISSVSTDANVAASKEIRNNDLRSSLSRNSFAPSQGSRDEYDTGSHGMSNLSSPSHVHESVTLAPLPQNPTIVGNRFTATSHHMHSTFDDKVLHRGNDADPLFLFAPATSHLRDVRKTSVVWDPEAGRYVSAPVTTTSEVRNRLLNPSSQTASTQNPRPILPAHDSSSGSSALRDPLPLHQAERRLTYTGDSIFYGGPLINIPTRDTPRSGRGLVRDVQDRLASTVHRDARIRRDSTSNQLPVFAPGGLGANSQTGSNIK",
            "Arabidopsis thaliana",
          ],
}

此后，您可以直接从哈希中查找序列和生物名称，而无需搜索。 sites 文件的读取循环就变成了这样。

while (<$sites_fh>) {

  my @chunks = split /\t/;
  next if $chunks[0] =~ /P3DB/;

  my ($pid, $acc, $psite) = $chunks[0, 1, 3];

  my $sequence = $proteins{$pid}[0];
  my $name     = $proteins{$pid}[1];

  # Use $sequence and $name somehow

}

请注意，我已将您的标识符更改为小写。如果您为包名称等全局变量保留大写字母，熟悉 Perl 的人会感谢您。

我还将您的文件句柄 AFILE 更改为词法 $sites_fh，因此您需要将 open 语句更改为类似

open my $sites_fh, '<', 'sites.pdb' or die $!;

我希望这会有所帮助。

【讨论】：