#!/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";
}
阅读(1322) | 评论(0) | 转发(0) |