Chinaunix首页 | 论坛 | 博客
  • 博客访问: 458733
  • 博文数量: 118
  • 博客积分: 4015
  • 博客等级: 上校
  • 技术积分: 1233
  • 用 户 组: 普通用户
  • 注册时间: 2010-11-24 22:11
文章分类

全部博文(118)

文章存档

2013年(5)

2011年(61)

2010年(52)

分类: Python/Ruby

2011-03-21 20:01:50

#!/usr/bin/perl
my $in1 = $ARGV[0];
my $in2 = $ARGV[1];
my $in3 = $ARGV[2];
die "NO INPUT 1!" if(!defined($in1));
die "NO INPUT 2!" if(!defined($in2));
die "NO INPUT 3!" if(!defined($in3));
open RH,$in1 or die "$!";
my @lines1 = ;
close RH or die "$!";

open RH,$in2 or die "$!";
my @lines2 = ; chomp @lines2;
close RH or die "$!";

open RH,$in3 or die "$!";
my @lines3 = ;
close RH or die "$!";

my $count = 0;
while($count < @lines1){
    my $minB = 100000000000;
    my $maxA ;
    my $gID = $lines2[$count];
    my $i = $count;
    for($i = $count; $i < @lines2 && $lines2[$i] == $gID; $i ++){
        my @fields = split(/\t/,$lines1[$i]);
    if($fields[7] < $minB){
        $minB = $fields[7];
    }
    $maxA = $fields[6];
    }
    my $fname = "";
    $lines1[$count] =~ /([0-9]+)-[0-9]+/;
    my $base = $1;
    my $st1 = $base + $maxA;
    my $st2 = $base + $minB;
    if($i - $count > 6){
        $fname = "can6_$st1-$st2.aln" ;   
    }else{
        $fname = "can_$st1-$st2.aln";
    }
    open WH,">$fname" or die "$!";


    my $text1 = "CLUSTAL W";
    $lines1[$count] =~ /(.*?)    $text1 = $text1."\tO_$1_($st1-$st2)";
    $lines1[$count] =~ /_(NC_\d+)_/;
    my $text2 = "$1_($st1-$st2)\t\t";
    my $text1E = "";
    my $countS = $count;
    my $z = 0;
    @fields = split(/\t/,$lines1[$count]);
    print "\t\tBEGIN LOOP1\n";
    for( $z = 0; $z < @lines3; $z ++){
        if($lines3[$z] =~ /Query= $fields[0]/){
        last;
    }
    }
    my $seqS = $fields[6];
    my $indexS = $maxA - $fields[6];
    my $seqlen = $minB - $maxA + 1;
    my $y = $z + 1;
L1:    while($lines3[$y] !~ /^Query=/){
          while($lines3[$y] !~ />$fields[1]/){
          $y ++;
      }
      $y ++;
      while($lines3[$y] !~ /^>/){
          if($lines3[$y] =~ /^Query:/){
              @arrs = split(/\s+/,$lines3[$y]);chomp @arrs;
          if($arrs[1] == $seqS){
              #Head matched
              if($arrs[3] == $fields[7]){
                  #Tail Matched
              $text2 = $text2.substr($arrs[2],$indexS,$seqlen)."\n";
              last L1;

              }else{
                  #Tailed Missed
              if($arrs[3] < $fields[7]){
                  #Remain Some
                  $text2 = $text2.substr($arrs[2],$indexS,$arrs[3] - $arrs[1] + 1 - $indexS);
                  $seqlen = $seqlen - ($arrs[3] - $arrs[1] + 1 - $indexS);
                  $indexS = 0;
                  $seqS = $arrs[3] + 1;
                  
              }else{
                  #No left
                  $text2 = $text2.substr($arrs[2],$indexS,$seqlen)."\n";
                  last L1;
              }
              }
          }else{
              #Head missed
              ;
          }

          }else{
              #Sbject:
          ;
          }
          $y ++;
      }
          $y ++;
    }
    print "\t\tLOOP1 done!\n";
    for($count = $count; $count < $i ; $count ++){
       @fields = split(/\t/,$lines1[$count]); chomp @fields;
       $fields[1] =~ /(.*?)       $text1E = $text1E."\tO_M".($count - $countS + 1)."_Identity_$fields[2]"."\tO_M".($count - $countS + 1)."_E_$fields[10]"."\tO_M".($count - $countS + 1)."_Score_$fields[11]";
       print "\t\tLine ",$count - $countS ,"doing...\n";
               my $st3 = $fields[8];
               my $st4 = $fields[9];
               $fields[1] =~ /(\d+)-\d+/;
               my $base2 = $1;
               if($st3 > $st4){
                   $st3 = $st3 - ($maxA - $fields[6]);
               $st3 = $st3 + $base2;
               $st4 = $st3 - ($minB - $maxA);
               }else{
                   $st3 = $st3 + ($maxA - $fields[6]);
               $st3 = $st3 + $base2;
               $st4 = $st3 + ($minB - $maxA);
               }
               $fields[1] =~ /_(NC_\d+)_/;
               $text2 = $text2."$1_($st3-$st4)\t\t";
               print "HHHH\n $text2";
       $text1 = $text1."\tM".($count - $countS + 1)."_$1($st3-$st4)";
       $seqS = $fields[8];
       $indexS = $maxA - $fields[6];
       $seqlen = $minB - $maxA + 1;
       #Begin Text2
        $y = $z + 1;
L2:    while($lines3[$y] !~ /^Query=/){
           while($lines3[$y] !~ />$fields[1]/){
           $y ++;
       }
       $y ++;
       while($lines3[$y] !~ /^>/){
           if($lines3[$y] =~ /^Sbjct:/){
               #Sbjct
           @arrs = split(/\s+/,$lines3[$y]);chomp @arrs;
           if($arrs[1] == $seqS){
               #Head Match
               if($arrs[3] == $fields[9]){
                   #Tail Match
               $text2 = $text2.substr($arrs[2],$indexS,$seqlen)."\n";
               last L2;
               }else{
                   #Tail Miss
               if($fields[8] < $fields[9]){
                   # -+
                   if($arrs[3] < $fields[9]){
                       # Remain
                   $text2 = $text2.substr($arrs[2],$indexS,$arrs[3] - $arrs[1] + 1 - $indexS);
                   $seqS = $arrs[3] + 1;
                   $seqlen = $seqlen - ($arrs[3] - $arrs[1] + 1 - $indexS);
                   $indexS = 0;
                   }else{
                       # NO Left
                   $text2 = $text2.substr($arrs[2],$indexS,$seqlen)."\n";
                   last L2;
                   }
               }else{
                   # +-
                   if($arrs[3] < $fields[9]){
                       #No left
                   $text2 = $text2.substr($arrs[2],$indexS,$seqlen)."\n";
                   last L2;
                   }else{
                       #Remain
                   $text2 = $text2.substr($arrs[2],$indexS,$arrs[1] - $indexS - $arrs[3] + 1);
                   $seqS = $arrs[3] - 1;
                   $seqlen = $seqlen - ($arrs[1] - $indexS - $arrs[3] + 1);
                   $indexS = 0
                   }
               }
               }

           }else{
               #Head Miss
               ;
           
           }

           }else{
               #Query
           ;
           }
           $y ++;
       }
      
        }

    }
    $text1 = $text1.$text1E;
    print WH $text1."\n\n\n".$text2."\n";
    close WH or die "$!";
    print "No.$count is done!\n";

}

阅读(1240) | 评论(0) | 转发(0) |
0

上一篇:trunc.pl

下一篇:Perl 实现简单的进度条

给主人留下些什么吧!~~