免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
最近访问板块 发新帖
楼主: little_joe

’求助!!!求助!!!文件匹配输出特定的内容‘ [复制链接]

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-06-19 17:34 |显示全部楼层
回复 40# little_joe
再试下,
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. sub collect_data{
  5.     my (%ha, %hFilt, $id);
  6.     open (my $fhA, '<', shift);
  7.     while (defined (local $_ = <$fhA>)){
  8.         if (/\A(\D\S+)/){
  9.             %hFilt = () if (!exists $ha{$1});
  10.             $id = $1;
  11.             next;
  12.         }
  13.         s/\s+\z//;
  14.         push (@{$ha{$id}}, $_) if (!$hFilt{$_}++);
  15.     }
  16.     close ($fhA);
  17.     \%ha;
  18. }

  19. sub compare{
  20.     my ($id, $rh, @aData) = (@_);
  21.     my %hStat;
  22.     foreach my $v (@{$rh->{$id}}){
  23.         foreach my $V (grep {$v >= $_->[0] && $v <= $_->[2]} @aData){
  24.             my ($cnt, $leng) = (0, $v - $V->[0] + 1);
  25.             if ((local $_ = substr ($V->[1], 0, $leng)) =~ /-/){
  26.                 $cnt = s/-//g;
  27.                 $cnt += length ($1) if (substr ($V->[1], $leng) =~ /\A(-+)/);
  28.             }
  29.             my @aChars = split (//, $V->[4]);
  30.             next if ($aChars[$v - $V->[0] + $cnt] eq '-');
  31.             $hStat{$aChars[$v - $V->[0] + $cnt]}++;
  32.         }
  33.     }
  34.     if (%hStat){
  35.         print "$id\t";
  36.         print "$_:$hStat{$_}\t" for keys %hStat;
  37.         print "\n";
  38.     }
  39. }

  40. my $rha = collect_data (shift);

  41. open (my $fh, '<', shift);
  42. my ($Flag, $Offset, $ID, @aData) = 0;
  43. while (defined (local $_ = <$fh>)){
  44.     next if (/\A\s*\z/ or !/\A(?:Query|Sbjct)/);
  45.     if (/\AQuery=\h+(\S+)/){
  46.         compare ($ID, $rha, splice (@aData)) if (@aData);
  47.         exists ($rha->{$1}) ? do {($ID, $Flag, $Offset) = ($1, 1, -1)} : ($Flag = 0);
  48.         next;
  49.     }
  50.     next if (!$Flag);
  51.     if (/\A(Query|Sbjct)\h+(\d+)\h+(\H+)\h+(\d+)/){
  52.         push (@{$aData[$1 eq "Query" ? ++$Offset : $Offset]}, $2, $3, $4);
  53.     }
  54. }
  55. close ($fh);
  56. compare ($ID, $rha, splice (@aData)) if ($Flag);
复制代码
perl abc.pl a b
---------------------------------
NP_415088.1-1   W:4
NP_415560.1-1   S:1
YP_025310.1-1   V:1     S:1
YP_588459.1-1   N:1     K:1
YP_002791252.1-1        S:1
NP_418697.1-1   W:7     V:7

论坛徽章:
0
发表于 2016-06-19 17:52 |显示全部楼层
非常感谢!! 回复 41# sunzhiguolu


   

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-06-19 18:24 |显示全部楼层
本帖最后由 sunzhiguolu 于 2016-06-21 13:48 编辑

回复 42# little_joe
用这个处理一下你的实际文件, 应该没啥问题了.
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;

  4. sub collect_data{
  5.     my (%ha, %hFilt, $id);
  6.     open (my $fhA, '<', shift);
  7.     while (defined (local $_ = <$fhA>)){
  8.         if (/\A([a-z]\S+)/i){
  9.             %hFilt = () if (!exists $ha{$1});
  10.             $id = $1;
  11.             next;
  12.         }
  13.         push (@{$ha{$id}}, $1) if (/\A\s*(\d+)/ and !$hFilt{$1}++);
  14.     }
  15.     close ($fhA);
  16.     \%ha;
  17. }

  18. sub compare{
  19.     my ($id, $rh, @aData) = (@_);
  20.     my %hStat;
  21.     foreach my $v (@{$rh->{$id}}){
  22.         foreach my $V (grep {$v >= $_->[0] && $v <= $_->[2]} @aData){
  23.             my ($cnt, $leng) = (0, $v - $V->[0] + 1);
  24.             if ((local $_ = substr ($V->[1], 0, $leng)) =~ /-/){
  25.                 $cnt = s/-//g;
  26.                 $cnt += length ($1) if (substr ($V->[1], $leng - 2, 1) eq '-' and substr ($V->[1], $leng) =~ /\A(-+)/);
  27.             }
  28.             my @aChars = split (//, $V->[4]);
  29.             next if ($aChars[$v - $V->[0] + $cnt] eq '-');
  30.             $hStat{$aChars[$v - $V->[0] + $cnt]}++;
  31.         }
  32.     }
  33.     if (%hStat){
  34.         print "$id\t";
  35.         print "$_:$hStat{$_}\t" for keys %hStat;
  36.         print "\n";
  37.     }
  38. }

  39. my $rha = collect_data (shift);

  40. open (my $fh, '<', shift);
  41. my ($Flag, $Offset, $ID, @aData) = 0;
  42. while (defined (local $_ = <$fh>)){
  43.     next if (/\A\s*\z/ or !/\A(?:Query|Sbjct)/);
  44.     if (/\AQuery=\h+(\S+)/){
  45.         compare ($ID, $rha, splice (@aData)) if (@aData);
  46.         exists ($rha->{$1}) ? do {($ID, $Flag, $Offset) = ($1, 1, -1)} : ($Flag = 0);
  47.         next;
  48.     }
  49.     next if (!$Flag);
  50.     if (/\A(Query|Sbjct)\h+(\d+)\h+(\H+)\h+(\d+)/){
  51.         push (@{$aData[$1 eq "Query" ? ++$Offset : $Offset]}, $2, $3, $4);
  52.     }
  53. }
  54. close ($fh);
  55. compare ($ID, $rha, splice (@aData)) if ($Flag);
复制代码

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-06-19 18:43 |显示全部楼层
本帖最后由 sunzhiguolu 于 2016-06-19 18:51 编辑

最后再问你一个问题, 对于这种情况 如何处理: (假设 A 中数值 78)
以下为对应的 B 文件的 Query, Sbjct 行:
Query  59   FSGDMLPD-GINRRWLTPAB---LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113

Query  59   FSGDMLPDGINRRWLTPAB-LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113

Query  59   FSGDMLPD-GINRRWLTPABT-LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113

此时, 对于 Query, Sbjct 字符应该落在哪里? (回帖时请使用上面这个示例进行说明, 这样我容易理解一点) 只是验证下我对于你先前说的那样内容是否已经理解, 代码逻辑是否存在纰漏.

论坛徽章:
0
发表于 2016-06-20 17:27 |显示全部楼层
Query  59   FSGDMLPD-GINRRWLTPAB---LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113

Query  59   FSGDMLPDGINRRWLTPAB-LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113

Query  59   FSGDMLPD-GINRRWLTPABT-LRQPVVCLNRSGNMVGQSVLPMAG-LSGKKKATKNSLF  116
              G       IN R  T A S++ +P+          G+++    G    KK+  KNSL
Sbjct  61   APGWWDKPPLINARVETAATSRMFKPLW-------QHGRAICFADGWFEWKKEGDKNSLI  113
如果文件A数值为78则对应的Query行和Sbjct行应该是红色的字母,您的想法是正确的,谢谢您的热情帮助及指教,谢谢您!回复 44# sunzhiguolu


   

论坛徽章:
0
发表于 2016-06-21 09:53 |显示全部楼层
stop_pipei_hb-1.word (1.32 MB, 下载次数: 1)

论坛徽章:
7
2015年迎新春徽章
日期:2015-03-04 10:01:442015年亚洲杯之科威特
日期:2015-04-17 17:27:042015亚冠之阿尔萨德
日期:2015-07-22 15:25:4015-16赛季CBA联赛之吉林
日期:2016-04-11 13:57:30操作系统版块每日发帖之星
日期:2016-06-24 06:20:00每日论坛发贴之星
日期:2016-06-24 06:20:0015-16赛季CBA联赛之辽宁
日期:2017-05-09 12:31:45
发表于 2016-06-21 09:58 |显示全部楼层
顶贴,理解能力弱

论坛徽章:
0
发表于 2016-06-21 10:28 |显示全部楼层
谢谢关注!回复 47# 251744647


   

论坛徽章:
307
程序设计版块每周发帖之星
日期:2016-04-08 00:41:33操作系统版块每日发帖之星
日期:2015-09-02 06:20:00每日论坛发贴之星
日期:2015-09-02 06:20:00程序设计版块每日发帖之星
日期:2015-09-04 06:20:00每日论坛发贴之星
日期:2015-09-04 06:20:00每周论坛发贴之星
日期:2015-09-06 22:22:00程序设计版块每日发帖之星
日期:2015-09-09 06:20:00程序设计版块每日发帖之星
日期:2015-09-19 06:20:00程序设计版块每日发帖之星
日期:2015-09-20 06:20:00每日论坛发贴之星
日期:2015-09-20 06:20:00程序设计版块每日发帖之星
日期:2015-09-22 06:20:00程序设计版块每日发帖之星
日期:2015-09-24 06:20:00
发表于 2016-06-21 13:20 |显示全部楼层
回复 46# little_joe
再用 43 楼的代码试下,

   

论坛徽章:
7
巳蛇
日期:2013-11-28 09:22:59天秤座
日期:2014-10-25 15:40:452015年辞旧岁徽章
日期:2015-03-03 16:54:152015年迎新春徽章
日期:2015-03-04 09:53:172015亚冠之德黑兰石油
日期:2015-07-15 08:46:452015亚冠之平阳省
日期:2015-11-08 16:27:53白银圣斗士
日期:2015-11-14 09:58:12
发表于 2016-06-30 08:51 |显示全部楼层
回复 49# sunzhiguolu


    他的文件是BLAST软件包中的BLASTP程序的输出,BioPerl模块中的Bio::SearchIO可以处理这种输出数据,可以直接取序列匹配信息,以及对应的序列。序列提取出来以后,再找对应的氨基酸,就方便多了!

  1. PS C:\Users\liuha> perldoc Bio::SearchIO
  2. NAME
  3.     Bio::SearchIO - Driver for parsing Sequence Database Searches (BLAST,
  4.     FASTA, ...)

  5. SYNOPSIS
  6.        use Bio::SearchIO;
  7.        # format can be 'fasta', 'blast', 'exonerate', ...
  8.        my $searchio = Bio::SearchIO->new( -format => 'blastxml',
  9.                                          -file   => 'blastout.xml' );
  10.        while ( my $result = $searchio->next_result() ) {
  11.            while( my $hit = $result->next_hit ) {
  12.             # process the Bio::Search::Hit::HitI object
  13.                while( my $hsp = $hit->next_hsp ) {
  14.                 # process the Bio::Search::HSP::HSPI object
  15.                }
  16.            }
  17.        }

  18. DESCRIPTION
  19.     This is a driver for instantiating a parser for report files from
  20.     sequence database searches. This object serves as a wrapper for the
  21.     format parsers in Bio::SearchIO::* - you should not need to ever use
  22.     those format parsers directly. (For people used to the SeqIO system it,
  23.     we are deliberately using the same pattern).

  24.     Once you get a SearchIO object, calling next_result() gives you back a
  25.     Bio::Search::Result::ResultI compliant object, which is an object that
  26.     represents one Blast/Fasta/HMMER whatever report.
复制代码

评分

参与人数 1信誉积分 +10 收起 理由
sunzhiguolu + 10 开眼界了, 很给力!!!

查看全部评分

您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

  

北京盛拓优讯信息技术有限公司. 版权所有 京ICP备16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122 niuxiaotong@pcpop.com 17352615567
未成年举报专区
中国互联网协会会员  联系我们:huangweiwei@itpub.net
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP