’求助！！！求助！！！文件匹配输出特定的内容‘

sunzhiguolu

回复 50# b114213903
谢谢大神指点, 如果您有时间的话, 也请多多回帖 向您学习.

   

b114213903

回复 51# sunzhiguolu


    入门级小菜鸟，所谓的生物信息学员，跟类似的数据打交道比较多，经常用BioPerl模块而已……

little_joe

不好意思，又来麻烦您，您之前帮忙写的脚本非常好用，受益匪浅，我想在问您一下，要是在这两个文件事例中加一点文件限制条件，其它的匹配模式都不变，这样是不是直接加一些条件语句就好了？
例如：在id匹配上之后，读取文件b时对Identities = 49/56 (88%)即序列相似性大于90%才进行序列匹配输出，否则直接跳出忽略，
这样您写的程序我尝试修改了一下，结果功夫不到家，只好又来请教，
这是您上次写的，麻烦您在帮忙给看看
O(∩_∩)O谢谢
#!/usr/bin/perl
use strict;
use warnings;


sub collect_data{
    my (%ha, %hFilt, $id);
    open (my $fhA, '<', shift);
    while (defined (local $_ = <$fhA>){
        if (/\A([a-z]\S+)/i){
            %hFilt = () if (!exists $ha{$1});
            $id = $1;
            next;
        }
        push (@{$ha{$id}}, $1) if (/\A\s*(\d+)/ and !$hFilt{$1}++);
    }
    close ($fhA);
    \%ha;
}


sub compare{
    my ($id, $rh, @aData) = (@_);
    my %hStat;
    foreach my $v (@{$rh->{$id}}){
        foreach my $V (grep {$v >= $_->[0] && $v <= $_->[2]} @aData){
            my ($cnt, $leng) = (0, $v - $V->[0] + 1);
            if ((local $_ = substr ($V->[1], 0, $leng)) =~ /-/){
                $cnt = s/-//g;
                $cnt += length ($1) if (substr ($V->[1], $leng - 2, 1) eq '-' and substr ($V->[1], $leng) =~ /\A(-+)/);
            }
            my @aChars = split (//, $V->[4]);
            next if ($aChars[$v - $V->[0] + $cnt] eq '-');
            $hStat{$aChars[$v - $V->[0] + $cnt]}++;
        }
    }
    if (%hStat){
        print "$id\t";
        print "$_hStat{$_}\t" for keys %hStat;
        print "\n";
    }
}


my $rha = collect_data (shift);


open (my $fh, '<', shift);
my ($Flag, $Offset, $ID, @aData) = 0;
while (defined (local $_ = <$fh>){
    next if (/\A\s*\z/ or !/\A(?uery|Sbjct)/);
    if (/\AQuery=\h+(\S+)/){
        compare ($ID, $rha, splice (@aData)) if (@aData);
        exists ($rha->{$1}) ? do {($ID, $Flag, $Offset) = ($1, 1, -1)} : ($Flag = 0);
        next;
    }
    next if (!$Flag);
    if (/\A(Query|Sbjct)\h+(\d+)\h+(\H+)\h+(\d+)/){
        push (@{$aData[$1 eq "Query" ? ++$Offset : $Offset]}, $2, $3, $4);
    }
}
close ($fh);
compare ($ID, $rha, splice (@aData)) if ($Flag);

文件A:
WP_000636925.1-1
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
WP_000567766.1-1
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        77
        79
        79
        79
        80
        80
        80
        80
        80
        80
WP_000636925.1-1
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
        63
        62
WP_001341795.1-1
        24
        24
        24
        24
        24
文件B:
BLASTP 2.2.29+


Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A.
Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J.
Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of
protein database search programs", Nucleic Acids Res. 25:3389-3402.


Reference for composition-based statistics: Alejandro A. Schaffer,
L. Aravind, Thomas L. Madden, Sergei Shavirin, John L. Spouge, Yuri
I. Wolf, Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with

Query= WP_001386572.1-1

Length=21


***** No hits found *****



Lambda      K        H        a         alpha
   0.314    0.134    0.424    0.792     4.96

Gapped
Lambda      K        H        a         alpha    sigma
   0.267   0.0410    0.140     1.90     42.6     43.6

Effective search space used: 1206911076


Query= WP_000636925.1-1

Length=105
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

lcl|CU928161.2_prot_CAR03545.1_2173  [gene=ECS88_2260] [protein=c...    104   2e-29


>lcl|CU928161.2_prot_CAR03545.1_2173 [gene=ECS88_2260] [protein=conserved hypothetical protein] [protein_id=CAR03545.1]
[location=2240570..2240740]
Length=56

Score =   104 bits (260),  Expect = 2e-29, Method: Compositional matrix adjust.
Identities = 49/56 (88%), Positives = 52/56 (93%), Gaps = 0/56 (0%)

Query  8   LLFTIMRSLYLKIRFPSDLRSINALNYLQVPSIDWSDSGYVQHFINVIEKMLTRKK  63
           +LFTIMRSLYLK RFPSDL+SINALNYLQVPSIDWSDSGY+QHFINVIEKM T  K
Sbjct  1   MLFTIMRSLYLKSRFPSDLKSINALNYLQVPSIDWSDSGYIQHFINVIEKMPTTNK  56



Lambda      K        H        a         alpha
   0.325    0.135    0.407    0.792     4.96

Gapped
Lambda      K        H        a         alpha    sigma
   0.267   0.0410    0.140     1.90     42.6     43.6

Effective search space used: 1406358592
Query= WP_000636925.1-1

Length=105
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

lcl|CU928161.2_prot_CAR03545.1_2173  [gene=ECS88_2260] [protein=c...    104   2e-29


>lcl|CU928161.2_prot_CAR03545.1_2173 [gene=ECS88_2260] [protein=conserved hypothetical protein] [protein_id=CAR03545.1]
[location=2240570..2240740]
Length=56

Score =   104 bits (260),  Expect = 2e-29, Method: Compositional matrix adjust.
Identities = 49/56 (88%), Positives = 52/56 (93%), Gaps = 0/56 (0%)

Query  8   LLFTIMRSLYLKIRFPSDLRSINALNYLQVPSIDWSDSGYVQHFINVIEKMLTRKK  63
           +LFTIMRSLYLK RFPSDL+SINALNYLQVPSIDWSDSGY+QHFINVIEKM T  K
Sbjct  1   MLFTIMRSLYLKSRFPSDLKSINALNYLQVPSIDWSDSGYIQHFINVIEKMPTTNK  56



Lambda      K        H        a         alpha
   0.325    0.135    0.407    0.792     4.96

Gapped
Lambda      K        H        a         alpha    sigma
   0.267   0.0410    0.140     1.90     42.6     43.6

Effective search space used: 1406358592
回复 43# sunzhiguolu

sunzhiguolu

回复 53# little_joe

对于这句话, 如何理解?

在id匹配上之后，读取文件b时对Identities = 49/56 (88%)即序列相似性大于90%才进行序列匹配输出，否则直接跳出忽略

little_joe

文件a和文件b的匹配方式不变，只是之前当文件a与文件b匹配成功后，直接读取序列输出想要的结果，现在这些条件都不变，只需要在匹配成功后对文件b中的Identities = 49/56 (88%)进行筛选，这里是88%
而设置的筛选条件是90%，也就是只有Identities = 49/56 (88%）后面括号的数值大于90%才进行下一步序列读取匹配
黄色的您标出来的数字的结果就是后面括号的数值，它代表两条序列的相似性，49/56就是说总共序列是56个字母，完全一样的是49个，
这样两条序列的相似性就为88%，
谢谢您的回复~~~
回复 54# sunzhiguolu

sunzhiguolu

回复 55# little_joe

它代表两条序列的相似性，49/56就是说总共序列是56个字母，完全一样的是49个

这个 49 是怎么得来的?

little_joe

>lcl|CU928161.2_prot_CAR03545.1_2173 [gene=ECS88_2260] [protein=conserved hypothetical protein] [protein_id=CAR03545.1]
[location=2240570..2240740]
Length=56

Score =   104 bits (260),  Expect = 2e-29, Method: Compositional matrix adjust.
Identities = 49/56 (88%), Positives = 52/56 (93%), Gaps = 0/56 (0%)

Query  8   LLFTIMRSLYLKIRFPSDLRSINALNYLQVPSIDWSDSGYVQHFINVIEKMLTRKK  63
           +LFTIMRSLYLK RFPSDL+SINALNYLQVPSIDWSDSGY+QHFINVIEKM T  K
Sbjct  1   MLFTIMRSLYLKSRFPSDLKSINALNYLQVPSIDWSDSGYIQHFINVIEKMPTTNK  56

49就是两条序列比较得出的完全相同的部分，总共是49个字母，它们每个相互对应，如上面用酸橙色标出的部分中的字母总共就是49个，它们与Query和Sbjct 的序列相对应这里由于格式问题没对太清楚，您可以看看文件。例如
Query  122  PGAAALNGDYIWDMILIPVSSLILS  146
                  PGAAALNGDYIWDMILIPVSSLILS
Sbjct   122  PGAAALNGDYIWDMILIPVSSLILS  146
这个序列的Identities = 25/25 (100%)，25就是它们在BLAST（生物学中常用的软件，用于分析基因序列的相似性）的过程中匹配到的序列长度。回复 56# sunzhiguolu

sunzhiguolu

49就是两条序列比较得出的完全相同的部分，总共是49个字母，它们每个相互对应

大师, 到底是那两条序列. 先前好像没有提到两条序列的事啊, 也许我忘记了 麻烦提个醒.

little_joe

是我没说清楚地原因，两条序列就是Query和Sbjct两条互相对比的序列，Query  122  PGAAALNGDYIWDMILIPVSSLILS  146
                  PGAAALNGDYIWDMILIPVSSLILS
Sbjct   122  PGAAALNGDYIWDMILIPVSSLILS  146
例如这个序列的Identities = 25/25 (100%)
而上面的红色字母就是Query和Sbjct匹配上的相似序列
也就是说经过比对Query和Sbjct总共匹配到25个氨基酸（就是字母）
而且这25个氨基酸完全一样则Identities = 25/25 (100%)
这条序列Identities = 49/56 (88%)
Query  8   LLFTIMRSLYLKIRFPSDLRSINALNYLQVPSIDWSDSGYVQHFINVIEKMLTRKK  63
               +LFTIMRSLYLK RFPSDL+SINALNYLQVPSIDWSDSGY+QHFINVIEKM T  K
Sbjct  1   MLFTIMRSLYLKSRFPSDLKSINALNYLQVPSIDWSDSGYIQHFINVIEKMPTTNK  56
是因为总共匹配到56个氨基酸，但是其中有7个是（+和空格），
+代表两条序列的氨基酸虽然不同但是两个氨基酸性质相近，空格是为了让更多的氨基酸匹配上而插入的，
这些都是生物分析软件BLAST按照一定生物信息学原则自动进行的回复 58# sunzhiguolu

sunzhiguolu

修改两处代码:

1. #!/usr/bin/perl
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. sub collect_data{
   
6.     my (%ha, %hFilt, $id);
   
7.     open (my $fhA, '<', shift);
   
8.     while (defined (local $_ = <$fhA>)){
   
9.         if (/\A([a-z]\S+)/i){
   
10.             %hFilt = () if (!exists $ha{$1});
    
11.             $id = $1;
    
12.             next;
    
13.         }
    
14.         push (@{$ha{$id}}, $1) if (/\A\s*(\d+)/ and !$hFilt{$1}++);
    
15.     }
    
16.     close ($fhA);
    
17.     \%ha;
    
18. }
    
19. 
    
20. sub compare{
    
21.     my ($id, $rh, @aData) = (@_);
    
22.     my %hStat;
    
23.     foreach my $v (@{$rh->{$id}}){
    
24.         foreach my $V (grep {$v >= $_->[0] && $v <= $_->[2]} @aData){
    
25.             my ($cnt, $leng) = (0, $v - $V->[0] + 1);
    
26.             if ((local $_ = substr ($V->[1], 0, $leng)) =~ /-/){
    
27.                 $cnt = s/-//g;
    
28.                 $cnt += length ($1) if (substr ($V->[1], $leng - 2, 1) eq '-' and substr ($V->[1], $leng) =~ /\A(-+)/);
    
29.             }
    
30.             my @aChars = split (//, $V->[4]);
    
31.             next if ($aChars[$v - $V->[0] + $cnt] eq '-');
    
32.             $hStat{$aChars[$v - $V->[0] + $cnt]}++;
    
33.         }
    
34.     }
    
35.     if (%hStat){
    
36.         print "$id\t";
    
37.         print "$_:$hStat{$_}\t" for keys %hStat;
    
38.         print "\n";
    
39.     }
    
40. }
    
41. 
    
42. my $rha = collect_data (shift);
    
43. 
    
44. open (my $fh, '<', shift);
    
45. my ($Flag, $Offset, $ID, @aData) = 0;
    
46. while (defined (local $_ = <$fh>)){
    
47.     next if (/\A\s*\z/ or !/\A(?:Query|Sbjct|Identities)/);    # 添加一条判断
    
48.     if (/\AQuery=\h+(\S+)/){
    
49.         compare ($ID, $rha, splice (@aData)) if (@aData);
    
50.         exists ($rha->{$1}) ? do {($ID, $Flag, $Offset) = ($1, 1, -1)} : ($Flag = 0);
    
51.         next;
    
52.     }
    
53.     $Flag = 0 if (/\bIdentities.+?(\d+)\%/ and $1 < 90);    # 如果百分比小于 90%, 则忽略后面的统计
    
54.     next if (!$Flag);
    
55.     if (/\A(Query|Sbjct)\h+(\d+)\h+(\H+)\h+(\d+)/){
    
56.         push (@{$aData[$1 eq "Query" ? ++$Offset : $Offset]}, $2, $3, $4);
    
57.     }
    
58. }
    
59. close ($fh);
    
60. compare ($ID, $rha, splice (@aData)) if ($Flag);
    

复制代码