Chinaunix

标题: perl新手求教...求教......文件匹配输出指定内容，非常感谢 [打印本页]
作者: little_joe    时间: 2016-08-17 16:52
标题: perl新手求教...求教......文件匹配输出指定内容，非常感谢
perl新手求教...求教......文件匹配输出指定内容，非常感谢

文件a：
WP_001198722.1;YP_026238.1;
37
126
WP_000075876.1;
58
204
WP_000072067.1;YP_026241.1;
42
50
WP_000077872.1;YP_025308.1;
128
201
文件b：
>lcl|NC_012947.1_prot_WP_001198722.1_11 [gene=ECBD_RS00055] [protein=/inference=EXISTENCE: similar to AA sequence:RefSeq:WP_001198703.1; catalyzes the formation of D-glyceraldehyde 3-phosphate an>] [protein_id=WP_001198722.1] [location=10464..11081]
MQWQTKLPLIAILRGITPDEALAHVGAVIDAGFDAVEIPLNSPQWEQSIPAIVDAYGDKALIGAGTVLKP
EQVDALARMGCQLIVTPNIHSEVIRRAVGYGMTVCPGCATATEAFTALEAGAQALKIFPSSAFGPQYIKA
LKAVLPSDIAVFAVGGVTPENLAQWIDAGCAGAGLGSDLYRAGQSVERTAQQAAAFVKAYREAVQ
>lcl|NC_012947.1_prot_WP_000072067.1_4 [gene=gyrB] [protein=/inference=EXISTENCE: similar to AA sequence:RefSeq:WP_001341774.1; negatively supercoils closed circular double-stranded DNA] [protein_id=WP_000072067.1] [location=3957..6371]
MSNSYDSSSIKVLKGLDAVRKRPGMYIGDTDDGTGLHHMVFEVVDNAIDEALAGHCKEIIVTIHADNSVS
VQDDGRGIPTGIHPEEGVSAAEVIMTVLHAGGKFDDNSYKVSGGLHGVGVSVVNALSQKLELVIQREGKI
HRQIYEHGVPQAPLAVTGETEKTGTMVRFWPSLETFTNVTEFEYEILAKRLRELSFLNSGVSIRLRDKRD
GKEDHFHYEGGIKAFVEYLNKNKTPIHPNIFYFSTEKDGIGVEVALQWNDGFQENIYCFTNNIPQRDGGT
HLAGFRAAMTRTLNAYMDKEGYSKKAKVSATGDDAREGLIAVVSVKVPDPKFSSQTKDKLVSSEVKSAVE
QQMNELLAEYLLENPTDAKIVVGKIIDAARAREAARRAREMTRRKGALDLAGLPGKLADCQERDPALSEL
YLVEGDSAGGSAKQGRNRKNQAILPLKGKILNVEKARFDKMLSSQEVATLITALGCGIGRDEYNPDKLRY
HSIIIMTDADVDGSHIRTLLLTFFYRQMPEIVERGHVYIAQPPLYKVKKGKQEQYIKDDEAMDQYQISIA
LDGATLHTNASAPALAGEALEKLVSEYNATQKMINRMERRYPKAMLKELIYQPTLTEADLSDEQTVTRWV
NALVSELNDKEQHGSQWKFDVHTNAEQNLFEPIVRVRTHGVDTDYPLDHEFITGGEYRRICTLGEKLRGL
LEEDAFIERGERRQPVASFEQALDWLVKESRRGLSIQRYKGLGEMNPEQLWETTMDPESRRMLRVTVKDA
IAADQLFTTLMGDAVEPRRAFIEENALKAANIDI
>lcl|NC_012947.1_prot_WP_000772941.1_8 [gene=ECBD_RS00040] [protein=/inference=EXISTENCE: similar to AA sequence:RefSeq:WP_001176565.1] [protein_id=WP_000772941.1] [location=complement(7982..863]
MKLNFKGFFKAAGLFPLALMLSGCISYALVSHTAKGSSGKYQSQSDTITGLSQAKDSNGTKGYVFVGESL
DYLITDGADDIVKMLNNPALNRHNIQVADDARFVLNAGKKKFTGTISLYYYWNNEEEKALATHYGFACGV
QHCTRSLENLKGTIHEKNKNMDYSKVMAFYHPFKVRFYEYYSPRGIPDGVSAALLPVTVTLDIITAPLQF，
LVVYAVNQ
如上，文件a中的WP_001198722.1匹配到文件b中的WP_001198722.1（均用红色注明）
则用文件a：WP_001198722.1下面的数字37,126查找文件b中的序列，序列的开头为青色标明处（该处序号记为1）
文件b每个以>lcl开头的为一条完整的注释序列，
图中橙色的部分为注释信息（需要用文件a的WP_001198722.1匹配文件b注释信息中的WP_001198722.1），
下面的大写字母为序列信息，即需要进行序号匹配的内容，例如WP_001198722.1下的37,126分别对应文件b中的E,K，

输出结果：
WP_001198722.1 E:1 K:1
WP_000072067.1 E:2
即对匹配到的序号输出的字母进行计数；
谢谢您的关注及帮助；

作者: sunzhiguolu    时间: 2016-08-17 19:48
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  5. my ($id, %hData, %hCnt);
  6. open (my $FHa, '<', './a');
  7. while (defined (local $_ = <$FHa>)){
  8.     if (/\A(WP_[^;]+)/){
  9.         $id = $1;
  10.         next;
  11.     }
  12.     chomp;
  13.     push (@{$hData{$id}}, $_);
  14. }
  15. close ($FHa);

  17. open (my $FHb, '<', './b');
  18. local $/ = ">lcl";
  19. while (defined (local $_ = <$FHb>)){
  20.     chomp;
  21.     next if (/\A\s*\z/);
  22.     my @aLines = split (/\n/);
  23.     if ($aLines[0] =~ /(WP_[\d.]+)/ and exists $hData{$1}){
  24.         my $context = join ("", @aLines[1 .. $#aLines]);
  25.         $hCnt{substr ($context, $_ - 1, 1)}++ for @{$hData{$1}};
  26.         my @aK2V = map {"$_:$hCnt{$_}"} keys %hCnt;
  27.         print "$1 @aK2V\n";
  28.         %hCnt = ();
  29.     }
  30. }
  31. close ($FHb);
复制代码
作者: little_joe    时间: 2016-08-18 07:55
非常感谢您的帮助，我将您写的代码命名为new1.txt，用事例运行之后出现这样的提示 ，没有看懂，希望您能帮忙解答一下，另外那个b文件实际上特别大的，这样会影响代码运行么？ 回复 2# sunzhiguolu


   
作者: sunzhiguolu    时间: 2016-08-18 11:39
你的测试环境是啥呀?
作者: sunzhiguolu    时间: 2016-08-18 12:56
再试下,
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  5. my ($id, %hData, %hCnt);
  6. open (my $FHa, '<', shift);
  7. while (defined (local $_ = <$FHa>)){
  8.     if (/\A(WP_[^;]+)/){
  9.         $id = $1;
  10.         next;
  11.     }
  12.     s/\A\s+|\s+\z//g;
  13.     push (@{$hData{$id}}, $_);
  14. }
  15. close ($FHa);

  17. open (my $FHb, '<', shift);
  18. local $/ = ">lcl";
  19. while (defined (local $_ = <$FHb>)){
  20.     chomp;
  21.     next if (!length);
  22.     my @aLines = split (/\r?\n/);
  23.     if ($aLines[0] =~ /(WP_[\d.]+)/ and exists $hData{$1}){
  24.         my $context = join ("", @aLines[1 .. $#aLines]);
  25.         $hCnt{substr ($context, $_ - 1, 1)}++ for @{$hData{$1}};
  26.         my @aK2V = map {"$_:$hCnt{$_}"} keys %hCnt;
  27.         %hCnt = ();
  28.         print "$1 @aK2V\n";
  29.     }
  30. }
  31. close ($FHb);
复制代码
perl abc.pl a.txt b.txt
---------------------------------------------------------------------------
WP_001198722.1 K:1 E:1
WP_000072067.1 E:2
作者: sunzhiguolu    时间: 2016-08-18 13:00
open (my $FHa, '<', './a');

我的上个代码将文件路径写死了, 改成 shift 就可以了!
作者: little_joe    时间: 2016-08-18 17:11
恩恩，谢谢您的热情帮助，有问题再问您，再次感谢 回复 6# sunzhiguolu


   
作者: little_joe    时间: 2016-08-18 17:12
应该是UNIX的回复 4# sunzhiguolu


   
作者: sunzhiguolu    时间: 2016-09-01 22:58
本帖最后由 sunzhiguolu 于 2016-09-01 23:52 编辑

还真是有点记不住题目的要求,
作者: sunzhiguolu    时间: 2016-09-01 23:02
本帖最后由 sunzhiguolu 于 2016-09-01 23:49 编辑

试下,
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  5. my ($id, %hData, %hCnt);
  6. open (my $FHr, '<', shift);
  7. while (local $_ = <$FHr>){
  8.     s/\s+\z//;
  9.     next if (!length);
  10.     if (/\A([a-z]+_[^;]+)/i){
  11.         $id = $1;
  12.         next;
  13.     }
  14.     push (@{$hData{$id}}, $_);
  15. }
  16. close ($FHr);

  18. local $/ = '>';
  19. open ($FHr, '<', shift);
  20. while (local $_ = <$FHr>){
  21.     chomp;
  22.     next if (!length);
  23.     my @aSeq = split (/-\d+;\r?\n/);
  24.     next if (!exists $hData{$aSeq[0]});
  25.     print $aSeq[0];
  26.     $hCnt{substr ($aSeq[-1], $_ - 2, 1)}++ for @{$hData{$aSeq[0]}};
  27.     print "\t$_:$hCnt{$_}" for keys %hCnt;
  28.     print "\n";
  29.     %hCnt = ();
  30. }
  31. close ($FHr);
复制代码
perl abc.pl a b
------------------------------------
NP_414552.1     *:1
NP_414557.1     R:1
NP_414562.1     A:1     M:12


作者: little_joe    时间: 2016-09-02 08:44
回复 10# sunzhiguolu
又麻烦您了，O(∩_∩)O谢谢，

作者: 华小飞_Perl    时间: 2016-09-20 21:33
写一个~
  1. #!/usr/bin/perl

  3. use warnings;
  4. use strict;

  6. open DATA_1, '<', 'a.txt' or die "Can't open: $!";
  7. open DATA_2, '<', 'b.faa' or die "Can't open: $!";

  9. my (%pro_seqID, $seqID);
  10. while (<DATA_1>) {
  11.         chomp;
  12.         if (/(WP.*?);/) {
  13.                 $seqID = $1;
  14.                 next;
  15.         }
  16.         push @{ $pro_seqID{$seqID} }, $_;
  17. }
  18. close DATA_1;

  20. local $/ = "\n>";
  21. while (<DATA_2>) {
  22.         chomp;
  23.         s/^>//;
  24.         if ( my @id = m/_(WP.*?\.\d)_/ ) {
  25.                 my $id = shift @id;
  26.                 s/^.*\]//;                       
  27.                 my $pro_seq = join "", split;                       
  28.                 if ( exists $pro_seqID{$id} ) {
  29.                         my %cal_str;
  30.                         my $part_str = substr( $pro_seq, ( $pro_seqID{$id}->[0] - 1 ),
  31.                                 ( $pro_seqID{$id}->[1] - $pro_seqID{$id}->[0] ) +1 );
  32.                         my @chr = $part_str =~ m/^(\w).*(\w)$/;
  33.                         $cal_str{$_}++ for @chr;
  34.                         print $id, " ", $_, ": ", $cal_str{$_}, "\n" for ( sort keys %cal_str );
  35.                 }
  36.         }
  37.         print "\n";
  38. }
  39. close DATA_2;
复制代码

处理fasta文件如果大同小异的话，直接写个子程序，以后方便自己调用就好~
作者: b114213903    时间: 2016-09-21 09:04
回复 13# 华小飞_Perl

处理fasta格式有现成的包  BioPerl中的Bio::SeqIO模块就可以了
作者: 华小飞_Perl    时间: 2016-09-21 13:49
回复 14# b114213903

听过Bioperl，有的话就直接拿来用了，谢谢提醒。



欢迎光临 Chinaunix (http://bbs.chinaunix.net/) Powered by Discuz! X3.2