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

little_joe

这些冗余数据可以去掉的，不好意思，一开始我没想到是不是应该去掉，还有每个文件A对应的ID匹配文件B后需要读取该ID下面的所有Query和Sbjct行，如果对应的Sbjct行不是字母是‘-'的话就忽略，麻烦您了回复 20# sunzhiguolu


   

sunzhiguolu

回复 21# little_joe
请看下 19 楼, 我的问题.

   

little_joe

第二行的内容不是Query行的衔接，它们是没有关系的，第二行其实是Query行和Sbjct行的关系描述，+表示两个氨基酸是同义氨基酸，这是生物学相关的问题，所以读取文件B时不用考虑第二行，只需要考虑以Query和Sbjct开头的行就行。
文件B中如您所标黄色的序号不会出现大的在前，小的在后的情况，因为氨基酸序列是连续的有顺序的。
因为每个Query行开头都会有一个数字序号，如果文件A中的数字小于文件BQuery第一个数字，那么便可以直接忽略，即如果文件A数字不在该区间内即可忽略，
如果文件A为35则直接忽略以2开头的Query行，只读取含有35的Query行回复 19# sunzhiguolu


   

sunzhiguolu

回复 23# little_joe
好的, 如有不明白 我再问.

   

little_joe

恩恩，您辛苦啦！谢谢！回复 24# sunzhiguolu


   

sunzhiguolu

回复 25# little_joe
试下,

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(\D\S+)/){
   
10.             %hFilt = () if (!exists $ha{$1});
    
11.             $id = $1;
    
12.             next;
    
13.         }
    
14.         s/\s+\z//;
    
15.         push (@{$ha{$id}}, $_) if (!$hFilt{$_}++);
    
16.     }
    
17.     close ($fhA);
    
18.     \%ha;
    
19. }
    
20. 
    
21. sub compare{
    
22.     my ($id, $rh, @aData) = (@_);
    
23.     my %hStat;
    
24.     foreach my $v (@{$rh->{$id}}){
    
25.         foreach my $V (grep {$v >= $_->[0] && $v <= $_->[2] and $v >= $_->[3] && $v <= $_->[-1]} @aData){
    
26.             my ($q1, $q2, $q3, $s1, $s2, $s3) = @$V;
    
27.             my @aChars = split //, ($s2 =~ s/-//gr);
    
28.             $hStat{$aChars[$v - $s1]}++;
    
29.         }
    
30.     }
    
31.     if (%hStat){
    
32.         print "$id\n";
    
33.         print "$_:$hStat{$_}\n" for keys %hStat;
    
34.     }
    
35. }
    
36. 
    
37. my $rha = collect_data (shift);
    
38. 
    
39. open (my $fh, '<', shift);
    
40. my ($Flag, $Offset, $ID, @aData) = 0;
    
41. while (defined (local $_ = <$fh>)){
    
42.     next if (/\A\s*\z/ or !/\A(?:Query|Sbjct)/);
    
43.     if (/\AQuery=\h+(\S+)/){
    
44.         compare ($ID, $rha, splice (@aData)) if (@aData);
    
45.         exists ($rha->{$1}) ? do {($ID, $Flag, $Offset) = ($1, 1, -1)} : ($Flag = 0);
    
46.         next;
    
47.     }
    
48.     next if (!$Flag);
    
49.     if (/\A(Query|Sbjct)\h+(\d+)\h+(\H+)\h+(\d+)/){
    
50.         $Offset++ if ($1 eq "Query");
    
51.         push (@{$aData[$Offset]}, $2, $3, $4);
    
52.     }
    
53. }
    
54. close ($fh);
    
55. compare ($ID, $rha, splice (@aData)) if ($Flag);
    

复制代码

perl abc.pl a b
----------------------------------------------------------------------------------------------
NP_414894.2-1
M:4
NP_415921.2-1
M:4
YP_025310.1-1
A:2
YP_026163.2-1
M:4
YP_002791252.1-1
C:1
NP_417338.3-1
M:4
NP_417516.3-1
M:4
NP_418692.2-1
M:4
NP_418697.1-1
L:1
H:1



   

stanley_tam

感觉我的结果不一样

1. #！perl
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. package Sequence;
   
6. 
   
7. sub new {
   
8.     my ($class, $id) = @_;
   
9. 
   
10.     my $self             = +{};
    
11.     $self->{'id'}        = $id;
    
12.     $self->{'members'}   = [];
    
13. 
    
14.     bless $self, $class;
    
15.     return $self;
    
16. }
    
17. 
    
18. sub add_member {
    
19.     my ($self, $member) = @_;
    
20.     push @{ $self->{'members'} }, $member;
    
21. }
    
22. 
    
23. sub get_members {
    
24.     my ($self) = @_;
    
25.     return @{ $self->{'members'} };
    
26. }
    
27. 
    
28. sub get_id {
    
29.     my ($self) = @_;
    
30.     return $self->{'id'};
    
31. }
    
32. 
    
33. 1;
    
34. 
    
35. package SequenceCollection;
    
36. 
    
37. sub new {
    
38.     my ($class)                = @_;
    
39.     my $self                   = +{};
    
40.     $self->{'sequences'}       = [];
    
41.     $self->{'_id_seq'}         = +{};
    
42.     $self->{'last_seen_id'}    = q{};
    
43.     $self->{'is_duplicate_id'} = 0;
    
44. 
    
45.     bless $self, $class;
    
46.     return $self;
    
47. }
    
48. 
    
49. sub process_line {
    
50.     my ($self, $line) = @_;
    
51. 
    
52.     $line =~ s{\s}{}gmix;
    
53.     return if not $line;
    
54. 
    
55.     if ($line =~ m{^ \d+ $}mix) {
    
56.         return if $self->{'is_duplicate_id'};
    
57.         my $last_seen_id = $self->{'last_seen_id'};
    
58.         my $sequence     = $self->{'_id_seq'}->{$last_seen_id};
    
59.         $sequence->add_member($line);
    
60.     }
    
61.     else {
    
62.         $self->{'last_seen_id'} = $line;
    
63. 
    
64. 
    
65.         if (exists $self->{'_id_seq'}->{$line}) {
    
66.             $self->{'is_duplicate_id'} = 1;
    
67.         }
    
68.         else {
    
69.             # new record
    
70.             $self->{'is_duplicate_id'} = 0;
    
71. 
    
72.             my $sequence = Sequence->new($line);
    
73.             $self->{'_id_seq'}->{$line} = $sequence;
    
74. 
    
75.             push @{$self->{'sequences'}}, $sequence;
    
76.         }
    
77.     }
    
78. }
    
79. 
    
80. sub get_sequences {
    
81.     my ($self) = @_;
    
82.     return @{ $self->{'sequences'} };
    
83. }
    
84. 
    
85. 1;
    
86. 
    
87. package Query;
    
88. 
    
89. sub new {
    
90.     my ($class, $id)                  = @_;
    
91.    
    
92.     my $self                          = +{};
    
93.     $self->{'id'}                     = $id;
    
94.     $self->{'query_list'}             = [];
    
95.     $self->{'subject_list'}           = [];
    
96.     $self->{'last_seen_start_number'} = undef;
    
97. 
    
98.     bless $self, $class;
    
99.     return $self;
    
100. }
     
101. 
     
102. sub populate_list {
     
103.     my ($self, $string_type, $string, $start_number) = @_;
     
104. 
     
105.     my @letters = split //, $string;
     
106. 
     
107.     if ($string_type eq 'Query') {
     
108.         my $query_list = $self->{'query_list'};
     
109.         $self->{'last_seen_start_number'} = $start_number;
     
110. 
     
111.         for my $letter (@letters){
     
112.             $query_list->[$start_number] = $letter;
     
113.             ++$start_number;
     
114.         }
     
115. 
     
116.         $self->{'query_list'} = $query_list;
     
117.     }
     
118.     elsif ($string_type eq 'Sbjct') {
     
119.         my $subject_list = $self->{'subject_list'};
     
120.         $start_number = $self->{'last_seen_start_number'};
     
121. 
     
122.         for my $letter (@letters){
     
123.             $subject_list->[$start_number] = $letter;
     
124.             ++$start_number;
     
125.         }
     
126. 
     
127.         $self->{'subject_list'} = $subject_list;
     
128.     }
     
129.     else {
     
130.         die "This shouldn't happen...$/";
     
131.     }
     
132. }
     
133. 
     
134. sub get_subject_character {
     
135.     my ($self, $number) = @_;
     
136.     my $subject_list    = $self->{'subject_list'};
     
137.     my $character       = $subject_list->[$number] // q{};
     
138.     return $character
     
139. }
     
140. 
     
141. 1;
     
142. 
     
143. package QueryCollection;
     
144. 
     
145. sub new {
     
146.     my ($class, $id)          = @_;
     
147.    
     
148.     my $self                  = +{};
     
149.     $self->{'queries'}        = [];
     
150.     $self->{'id_query'}       = +{};
     
151.     $self->{'last_seen_id'}   = undef;
     
152. 
     
153.     bless $self, $class;
     
154.     return $self;
     
155. }
     
156. 
     
157. sub process_line {
     
158.     my ($self, $line) = @_;
     
159.     $line             =~ s{^\s+ | \s+$}{}gmix;
     
160. 
     
161.     return if not $line;
     
162. 
     
163.     if ($line =~ m{Query=}) {
     
164.         my ($id)                = $line =~ m{^Query= \s* (.*)$}mx;
     
165.         $self->{'last_seen_id'} = $id;
     
166. 
     
167.         my $query = Query->new($id);
     
168. 
     
169.         push @{ $self->{'queries'} }, $query;
     
170.         $self->{'id_query'}->{$id}  = $query;
     
171.     }
     
172.     elsif ($line =~ m{Query \s+ \d+ \s+}mix){
     
173.         my $last_seen_id            = $self->{'last_seen_id'};
     
174.         my $query                   = $self->{'id_query'}->{$last_seen_id};
     
175.         my ($start_number, $string) = $line =~ m{^Query \s+ (\d+) \s+ (\S+) \s+}mx;
     
176.         my $string_type             = 'Query';
     
177.         $query->populate_list($string_type, $string, $start_number);
     
178.     }
     
179.     elsif ($line =~ m{Sbjct \s+ \d+ \s+}mix){
     
180.         my $last_seen_id            = $self->{'last_seen_id'};
     
181.         my $query                   = $self->{'id_query'}->{$last_seen_id};
     
182.         my ($start_number, $string) = $line =~ m{^Sbjct \s+ (\d+) \s+ (\S+) \s+}mx;
     
183.         my $string_type             = 'Sbjct';
     
184.         $query->populate_list($string_type, $string, $start_number);
     
185.     }
     
186. }
     
187. 
     
188. sub get_query_by_id {
     
189.     my ($self, $id) = @_;
     
190.     my $query       = $self->{'id_query'}->{$id} || q{};
     
191.     return $query;
     
192. }
     
193. 
     
194. sub get_subject_character {
     
195.     my ($self, $id, $number) = @_;
     
196.     my $character            = q{};
     
197.    
     
198.     my $query                = $self->get_query_by_id($id);
     
199. 
     
200.     if ($query) {
     
201.         $character = $query->get_subject_character($number);
     
202.     }
     
203. 
     
204.     return $character;
     
205. }
     
206. 
     
207. 1;
     
208. 
     
209. package main;
     
210. 
     
211. 
     
212. sub main {
     
213.     my $sequence_collection = SequenceCollection->new();
     
214.     open my $a_fh, '<', 'A.word';
     
215.     while (defined(my $line = readline $a_fh)) {
     
216.         $sequence_collection->process_line($line);
     
217.     }
     
218.     close $a_fh;
     
219. 
     
220. 
     
221.     my $query_collection = QueryCollection->new();
     
222.     open my $b_fh, '<', 'B.word';
     
223.     while (defined(my $line = readline $b_fh)) {
     
224.         $query_collection->process_line($line);
     
225.     }
     
226.     close $b_fh;
     
227. 
     
228. 
     
229.     for my $sequence ( $sequence_collection->get_sequences() ){
     
230.         my @members = $sequence->get_members();
     
231.         my $id          = $sequence->get_id();
     
232. 
     
233.         my %count = ();
     
234.         for my $number ( @members ){
     
235.             my $character = $query_collection->get_subject_character($id, $number);
     
236.             ++$count{$character};
     
237.         }
     
238. 
     
239.         print "$/id => [$id]$/";
     
240.         print "members => [@members]$/";
     
241.         for my $character (keys %count){
     
242.             my $number = $count{$character};
     
243.             print "character => [$character]$/";
     
244.             print "number => [$number]$/";
     
245.         }
     
246. 
     
247.     }
     
248. 
     
249. }
     
250. 
     
251. main();
     
252. 
     
253. 
     
254. __END__

复制代码

输出：

1. id => [NP_414894.2-1]
   
2. members => [77]
   
3. character => [-]
   
4. number => [1]
   
5. 
   
6. id => [NP_415088.1-1]
   
7. members => [134]
   
8. character => [W]
   
9. number => [1]
   
10. 
    
11. id => [NP_415560.1-1]
    
12. members => [137]
    
13. character => [S]
    
14. number => [1]
    
15. 
    
16. id => [NP_415921.2-1]
    
17. members => [77]
    
18. character => [-]
    
19. number => [1]
    
20. 
    
21. id => [YP_025310.1-1]
    
22. members => [73 78]
    
23. character => [T]
    
24. number => [1]
    
25. character => [V]
    
26. number => [1]
    
27. 
    
28. id => [YP_026163.2-1]
    
29. members => [77]
    
30. character => [-]
    
31. number => [1]
    
32. 
    
33. id => [YP_588459.1-1]
    
34. members => [62 63]
    
35. character => [N]
    
36. number => [1]
    
37. character => [K]
    
38. number => [1]
    
39. 
    
40. id => [YP_002791252.1-1]
    
41. members => [26]
    
42. character => [S]
    
43. number => [1]
    
44. 
    
45. id => [NP_417338.3-1]
    
46. members => [77]
    
47. character => [-]
    
48. number => [1]
    
49. 
    
50. id => [NP_417516.3-1]
    
51. members => [77]
    
52. character => [-]
    
53. number => [1]
    
54. 
    
55. id => [NP_418692.2-1]
    
56. members => [77]
    
57. character => [-]
    
58. number => [1]
    
59. 
    
60. id => [NP_418697.1-1]
    
61. members => [81 85]
    
62. character => [W]
    
63. number => [1]
    
64. character => [E]
    
65. number => [1]
    

复制代码

sunzhiguolu

stanley_tam 发表于 2016-06-18 16:34
感觉我的结果不一样输出：

差异之处, 比如:
A 文件:
NP_415088.1-1
134

对应 B 文件 Query= NP_415088.1-1
...
Query  134  W  134
            W
Sbjct  61   W  61
...
Query  134  W  134
            W
Sbjct  61   W  61
...
Query  134  W  134
            W
Sbjct  61   W  61

A 与 B 中的对应数值 134 在 B 中的 Query 134 W 134 中满足匹配要求, 但是其下的 Sbjct 61 W 61 行的 却不能满足 A 中数值 134 的匹配范围.
在这里我采用的是 双过滤 即 A 中的 134 必须在 B 中 Query, Sbjct 行中都满足要求才认为满足要求,  这一点可能需要楼主说明一下. (这个地方我的确拿不准)

little_joe

谢谢各位的关注和帮助，是我之前没说明白的原因，文件A的数值134其实是不用满足Sbjct的条件的，这里134的作用是找到对应Query行的对应字母，然后通过该字母（Query行）找到Sbjct行中Query行字母对应的字母（也就是其下方对应的字母），输出的是Sbjct行的字母，并且对这个输出的字母进行计数，文件中当Query行中出现“-”时不会影响序号区间例如：Query：1 ASDFGHJKL 9和Query：1  ASD-FGHJKL 9它们两个都是包含9个元素，“-”不会影响字母的序号，而当Query行是字母对应的Sbjct行是“-”时则不用输出例如：Query：1 ASDFGHJKL 9
                                                        Sbjct：1 ASD-FGHJK 8假设文件A数值为4，则找到了Query行的F此时Sbjct行对应“-”则可以不用输出（或者为了简单期间输出计数也可）这个应该是文件的输出结果（例子中较少，所以手动验证了一下）
NP_415088.1-1          W:4       
NP_415560.1-1                S:1       
YP_025310.1-1                S:1        V:1
NP_418697.1-1                W:7        V:7
YP_588459.1-1                N:1        K:1
YP_002791252.1-1        S:1       
回复 28# sunzhiguolu


   

little_joe

谢谢您，我之前有可能没说清楚，文件中当Query行中出现“-”时不会影响序号区间例如：Query：1 ASDFGHJKL 9和Query：1  ASD-FGHJKL 9它们两个都是包含9个元素，“-”不会影响字母的序号例如：Query：1  ASDF-GHJKL  9
                                            Sbjct：6  ASDFGIHLUL  16这时候Query行出现了“-”要是文件A数值为6时，则找到的Query行字母应该为H，对应的输出是Sbjct行的H,回复 27# stanley_tam