| 1 |
gpertea |
23 |
#!/usr/bin/perl |
| 2 |
|
|
|
| 3 |
|
|
#some filtering can be performed on the output |
| 4 |
|
|
# |
| 5 |
|
|
# this should be rewritten in C/C++ to be faster |
| 6 |
|
|
|
| 7 |
|
|
use strict; |
| 8 |
|
|
|
| 9 |
|
|
my $usage = q{ |
| 10 |
|
|
blast2btab [HITLEN=nn[%]] [OVHANG=nn[%]] [PID=nn] [PSIM=nn] |
| 11 |
|
|
[SCORE=nnnn | E=nnnnn] |
| 12 |
|
|
|
| 13 |
|
|
Format and filter wu-blast2 output coming at stdin. The output consists of |
| 14 |
|
|
the following tab delimited fields: |
| 15 |
|
|
0 1 2 3 4 5 |
| 16 |
|
|
q_name, date, q_len, method, db_file, hit_name, |
| 17 |
|
|
6 7 8 9 10 11 12 |
| 18 |
|
|
q_n5, q_n3, hit_n5, hit_n3, %id, %sim, score, |
| 19 |
|
|
13 14 15 16 |
| 20 |
|
|
file_pos1, file_pos2, hit_description, Frame, |
| 21 |
|
|
17 18 19 [ 20 21 ] |
| 22 |
|
|
strand, hit_len, E-value [ hit-coverage HSPs ] |
| 23 |
|
|
|
| 24 |
|
|
Filtering: |
| 25 |
|
|
PID for minimum percent identity |
| 26 |
|
|
PSIM for minimum percent similarity |
| 27 |
|
|
SCORE for minimum alignment score |
| 28 |
|
|
E for maximum E value (Blast E parameter) |
| 29 |
|
|
HITLEN for minimum overlap length |
| 30 |
|
|
OVHANG for maximum overhang length (mismatched overhangs) |
| 31 |
|
|
|
| 32 |
|
|
If '%' is specified for HITLEN or OVHANG values, the actual cutoff |
| 33 |
|
|
lengths will be computed as percentages of the query sequence length; |
| 34 |
|
|
|
| 35 |
|
|
If -topcomboX option is used, each line will represent a group, with the individual |
| 36 |
|
|
segment hits in the 21th field, and the actual hit coverage length as the 20th field |
| 37 |
|
|
}; |
| 38 |
|
|
|
| 39 |
|
|
die($usage) if ($ARGV[0]=~/\-h|\-?\-help/i); |
| 40 |
|
|
my $_bignum=9999999; |
| 41 |
|
|
my %a=('HITLEN'=>0, |
| 42 |
|
|
'OVHANG'=>$_bignum, |
| 43 |
|
|
'PID'=>0, |
| 44 |
|
|
'PSIM'=>0, |
| 45 |
|
|
'SCORE'=>0, |
| 46 |
|
|
'E'=>$_bignum, |
| 47 |
|
|
'DEFLINE'=>'none', |
| 48 |
|
|
'PARSE'=>0 |
| 49 |
|
|
); |
| 50 |
|
|
# try to parse filtering options, if any; |
| 51 |
|
|
my ($is_lenperc, $is_ovperc)=(0,0); |
| 52 |
|
|
my $is_filtered=0; |
| 53 |
|
|
my $bytes_read=0; |
| 54 |
|
|
foreach (@ARGV) { |
| 55 |
|
|
my @param=split(/=/); |
| 56 |
|
|
$is_filtered=1 if ($param[0] ne 'DEFLINE'); |
| 57 |
|
|
die $usage if (@param<2 || not exists($a{$param[0]})); |
| 58 |
|
|
if ($param[0] eq 'HITLEN') { |
| 59 |
|
|
my $p=rindex($param[1],'%'); |
| 60 |
|
|
if ($p>0) { |
| 61 |
|
|
$is_lenperc=1; |
| 62 |
|
|
$param[1]=substr($param[1],0,$p); |
| 63 |
|
|
} |
| 64 |
|
|
} |
| 65 |
|
|
elsif ($param[0] eq 'OVHANG') { |
| 66 |
|
|
my $p=rindex($param[1],'%'); |
| 67 |
|
|
if ($p>=0) { |
| 68 |
|
|
$is_ovperc=1; |
| 69 |
|
|
$param[1]=substr($param[1],0,$p); |
| 70 |
|
|
} |
| 71 |
|
|
} |
| 72 |
|
|
$a{$param[0]}=$param[1]; |
| 73 |
|
|
} |
| 74 |
|
|
my $defln=lc($a{'DEFLINE'}); |
| 75 |
|
|
if ($defln eq 'none') { |
| 76 |
|
|
$defln=0; |
| 77 |
|
|
} |
| 78 |
|
|
elsif ($defln eq 'q') { |
| 79 |
|
|
$defln=1; |
| 80 |
|
|
} |
| 81 |
|
|
elsif ($defln eq 'h') { |
| 82 |
|
|
$defln=2; |
| 83 |
|
|
} |
| 84 |
|
|
elsif ($defln eq 'b') { |
| 85 |
|
|
$defln=3; |
| 86 |
|
|
} |
| 87 |
|
|
else { |
| 88 |
|
|
die "$usage. Error at DEFLINE option."; |
| 89 |
|
|
} |
| 90 |
|
|
|
| 91 |
|
|
$/ = "\n"; |
| 92 |
|
|
#$\ = "\n"; |
| 93 |
|
|
#$, = "\t"; |
| 94 |
|
|
my ($minhlen, $maxoh, $doparse, $minpid, $minpsim, $minscore, $maxE)=( |
| 95 |
|
|
$a{'HITLEN'}, |
| 96 |
|
|
$a{'OVHANG'}, |
| 97 |
|
|
$a{'PARSE'}, |
| 98 |
|
|
$a{'PID'}, |
| 99 |
|
|
$a{'PSIM'}, |
| 100 |
|
|
$a{'SCORE'}, |
| 101 |
|
|
$a{'E'}); |
| 102 |
|
|
$maxE='1'.$maxE if lc(substr($maxE,0,1)) eq 'e'; |
| 103 |
|
|
|
| 104 |
|
|
my $writeline= ($is_filtered) ? \&write_flt : \&write_noflt; |
| 105 |
|
|
my $curgrp; #current group # |
| 106 |
|
|
my @grpdata; #when groups are defined, this accumulates HSP data |
| 107 |
|
|
# list of [$q5, $q3, $h5, $h3, $pid, $psim, $score, $frame] |
| 108 |
|
|
#my ($grphspcount, $grplen, $grpids, $grpsims); |
| 109 |
|
|
#segments expected, segments collected already, |
| 110 |
|
|
# total bases, total identities, total positives |
| 111 |
|
|
|
| 112 |
|
|
#only one file is accepted here... |
| 113 |
|
|
my ($query, $qlen, $q5, $q3, $hit, $hitlen, $h5, $h3, $ids, $sims, $bcount, $pid, $psim, |
| 114 |
|
|
$score, $evalue, $strand, $frame, $hdefline, $searchdb, $fpos, $fend, $qdefline); |
| 115 |
|
|
my $flush=0; |
| 116 |
|
|
$fpos=-1; |
| 117 |
|
|
$fend=-1; |
| 118 |
|
|
my $method='prog'; |
| 119 |
|
|
$searchdb='n/a'; |
| 120 |
|
|
my @lcdate=localtime(); |
| 121 |
|
|
my $date=($lcdate[4]+1).'/'.$lcdate[3].'/'.(1900+$lcdate[5]); |
| 122 |
|
|
$frame='*'; |
| 123 |
|
|
BIGLOOP: |
| 124 |
|
|
while (<STDIN>) { |
| 125 |
|
|
die "Line $. doesn't appear to be a blast output:\n$_\n" unless (/^(T?BLAST[NXP]?)/); |
| 126 |
|
|
NEXT_QRYBLOCK: |
| 127 |
|
|
$bytes_read+=length($_); |
| 128 |
|
|
next if /^\s*$/; |
| 129 |
|
|
|
| 130 |
|
|
goto QRYFOUND if m/^Query=\s*(\S+)\s*(.*)/; |
| 131 |
|
|
|
| 132 |
|
|
$method=lc($1); |
| 133 |
|
|
do { |
| 134 |
|
|
$_=<STDIN>; $bytes_read+=length($_); |
| 135 |
|
|
|
| 136 |
|
|
last BIGLOOP unless $_; |
| 137 |
|
|
} until m/^Query=\s*(\S+)\s*(.*)/; |
| 138 |
|
|
QRYFOUND: |
| 139 |
|
|
&write_group() if $query && $curgrp; |
| 140 |
|
|
$query=$1; |
| 141 |
|
|
$qdefline=''; |
| 142 |
|
|
if ($defln==1 || $defln==3) { |
| 143 |
|
|
$qdefline=$2; |
| 144 |
|
|
#print STDERR "qdefline set to '$qdefline'\n"; |
| 145 |
|
|
$qdefline =~ tr/\n//d; |
| 146 |
|
|
} |
| 147 |
|
|
my $linepos=$.; |
| 148 |
|
|
do { |
| 149 |
|
|
$_=<STDIN>; $bytes_read+=length($_); |
| 150 |
|
|
last BIGLOOP unless $_; |
| 151 |
|
|
&err_parse("Cannot find 'Length' line for Query $query") |
| 152 |
|
|
if ($.-$linepos) > 10000; |
| 153 |
|
|
} until /\s+\(([\d\,]+) letters/; |
| 154 |
|
|
|
| 155 |
|
|
$qlen=$1; |
| 156 |
|
|
$qlen =~ tr/,//d; |
| 157 |
|
|
# parse all matches for this $query: |
| 158 |
|
|
$_=<STDIN>; $bytes_read+=length($_); |
| 159 |
|
|
my $was_param; |
| 160 |
|
|
while ($_ && !(($was_param=m/^Parameters:/) || m/^EXIT CODE/)) { |
| 161 |
|
|
chomp; |
| 162 |
|
|
if (m/^Database:\s+(\S+)/) { |
| 163 |
|
|
$searchdb=$1; |
| 164 |
|
|
($searchdb)=($searchdb =~ m/([^\/]+)$/); |
| 165 |
|
|
} |
| 166 |
|
|
# Hit section: |
| 167 |
|
|
if (/^>(\S+)\s*(.*)/) { |
| 168 |
|
|
#HSPs here |
| 169 |
|
|
$hit=$1; |
| 170 |
|
|
#print STDERR "hitline is: $_\n"; |
| 171 |
|
|
#if ($defln>=2) { |
| 172 |
|
|
$hdefline=$2; |
| 173 |
|
|
#print STDERR "hdefline set to '$hdefline'\n"; |
| 174 |
|
|
# } |
| 175 |
|
|
NEXT_SUBJ: |
| 176 |
|
|
my $savepos; |
| 177 |
|
|
#&err_parse("Did you use '-noseqs' option for this blast search?") |
| 178 |
|
|
# if ($hit && $score && !$h3); |
| 179 |
|
|
$savepos=$.; |
| 180 |
|
|
while (<STDIN>) { |
| 181 |
|
|
$bytes_read+=length($_); |
| 182 |
|
|
&err_parse("Cannot find 'Length' line for hit $hit") |
| 183 |
|
|
if ($.-$savepos) > 10000; |
| 184 |
|
|
if (/^\s+Length\s+=\s+([\d\,]+)/) { |
| 185 |
|
|
$hitlen=$1; $hitlen =~ tr/,//d; |
| 186 |
|
|
last; |
| 187 |
|
|
} |
| 188 |
|
|
chomp; |
| 189 |
|
|
s/^\s+//; |
| 190 |
|
|
$hdefline.=' '.$_ if $_; |
| 191 |
|
|
} |
| 192 |
|
|
$h3=0;$score=0;$strand='+'; |
| 193 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 194 |
|
|
$savepos=$.+80; |
| 195 |
|
|
my $withseqs=0; |
| 196 |
|
|
while ($_) { |
| 197 |
|
|
#search of all HSPs until Subject match or EOF is encountered |
| 198 |
|
|
#print "&&&&&& $query: found hit $hit, score=$score\n"; |
| 199 |
|
|
if (/^>(\S+)\s*(.*)/) { |
| 200 |
|
|
#print "prevscore= $score, prevhit=$hit, prevcoords=$h3\n"; |
| 201 |
|
|
#&err_parse("-- coordinates not found for hit $hit. Please be sure blastn is given -noseqs option!") |
| 202 |
|
|
# if ($hit && $score && !$h3); |
| 203 |
|
|
&write_group() if $curgrp; |
| 204 |
|
|
chomp; |
| 205 |
|
|
$hit=$1;$hdefline=$2; |
| 206 |
|
|
goto NEXT_SUBJ; |
| 207 |
|
|
} |
| 208 |
|
|
goto SKIP_EMPTYLN if ($_ eq "\n"); |
| 209 |
|
|
if (($was_param=m/^Parameters:/) || m/^EXIT CODE/) { |
| 210 |
|
|
&write_group() if $curgrp; |
| 211 |
|
|
goto NEXT_QUERY; |
| 212 |
|
|
} |
| 213 |
|
|
if (/^\s*Score[\s=]+(\d+)/) { |
| 214 |
|
|
$score=$1; |
| 215 |
|
|
$fpos=$bytes_read-length($_); |
| 216 |
|
|
($evalue)=(m/\,\s*Expect\s*=\s*([\de\.\-]+)/); |
| 217 |
|
|
#print("got e-value $evalue for hit $hit \n"); |
| 218 |
|
|
my ($group)=(m/Group\s+=\s+(\d+)/); |
| 219 |
|
|
&write_group() if ($curgrp && $curgrp!=$group); |
| 220 |
|
|
$curgrp=$group; |
| 221 |
|
|
$savepos=$.+80; |
| 222 |
|
|
} |
| 223 |
|
|
elsif (/^\s*Identities\s+=\s+(\d+)\/(\d+)\s+\(([\d\.]+)\%\)/) { |
| 224 |
|
|
($ids, $bcount, $pid)=($1,$2,$3); |
| 225 |
|
|
$psim=$pid; |
| 226 |
|
|
if (/Positives\s+=\s+(\d+)\/\d+\s+\(([\d\.]+)\%\)/) { |
| 227 |
|
|
($sims, $psim)=($1, $2) |
| 228 |
|
|
} |
| 229 |
|
|
$strand=/Minus/?'Minus':'Plus'; |
| 230 |
|
|
if (/Frame\s+\=\s+([\+\-\d]+)/) { |
| 231 |
|
|
$frame=$1; |
| 232 |
|
|
$strand=(substr($frame,0,1) eq '-') ? 'Minus':'Plus'; |
| 233 |
|
|
} |
| 234 |
|
|
} |
| 235 |
|
|
elsif (m/^Query:\s+/) { #parse query coordinates |
| 236 |
|
|
if (m/^Query:\s*(\d+)[\s]+\-\s+(\d+)/) { #noseqs case |
| 237 |
|
|
($q5,$q3)=($1,$2); |
| 238 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 239 |
|
|
if (m/^Sbjct:\s*(\d+)\s+\-\s+(\d+)/) { #noseqs case |
| 240 |
|
|
($h5,$h3)=($1,$2); |
| 241 |
|
|
$fend=$bytes_read; |
| 242 |
|
|
if ($curgrp) { |
| 243 |
|
|
push(@grpdata, [$q5, $q3, $h5, $h3, $pid, $psim, $score, $frame, $ids, $sims, $bcount]); |
| 244 |
|
|
#$grphspcount++; |
| 245 |
|
|
#$grpids+=$ids; |
| 246 |
|
|
#$grpsims+=$sims; |
| 247 |
|
|
#$grplen+=$bcount; |
| 248 |
|
|
} |
| 249 |
|
|
else { |
| 250 |
|
|
&$writeline; |
| 251 |
|
|
} |
| 252 |
|
|
} |
| 253 |
|
|
else { |
| 254 |
|
|
&err_parse("Cannot find Sbjct: line for [noseqs] output (query=$query, subj=$hit)"); |
| 255 |
|
|
} |
| 256 |
|
|
} |
| 257 |
|
|
else { #read/parse the whole alignment to get the start-end coordinates |
| 258 |
|
|
#for both query and sequence |
| 259 |
|
|
if (m/^Query:\s+(\d+)\s+[A-Z,\-,\*]+\s+(\d+)$/i) { |
| 260 |
|
|
($q5,$q3)=($1,$2); |
| 261 |
|
|
} |
| 262 |
|
|
else { |
| 263 |
|
|
&err_parse("Invalid start of Query: alignment line (query=$query, subj=$hit)"); |
| 264 |
|
|
} |
| 265 |
|
|
$withseqs=1; |
| 266 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 267 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 268 |
|
|
if (m/^Sbjct:\s+(\d+)\s+[A-Z,\-,\*]+\s+(\d+)$/i) { |
| 269 |
|
|
($h5, $h3) = ($1, $2); |
| 270 |
|
|
} |
| 271 |
|
|
else { |
| 272 |
|
|
&err_parse("Cannot find start Sbjct: line for alignment (query=$query, subj=$hit)"); |
| 273 |
|
|
} |
| 274 |
|
|
$_=<STDIN>; $bytes_read+=length($_);#skip empty line |
| 275 |
|
|
#$_=<STDIN>; |
| 276 |
|
|
#print STDERR "Alignment continuation is:\n$_\n\n"; |
| 277 |
|
|
#exit(0); |
| 278 |
|
|
while (defined($_=<STDIN>) && m/^Query:\s+\d+\s+[A-Z,\-,\*]+\s+(\d+)$/i) { |
| 279 |
|
|
$bytes_read+=length($_); |
| 280 |
|
|
#print STDERR "found alignment continuation\n"; |
| 281 |
|
|
$q3=$1; |
| 282 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 283 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 284 |
|
|
if (m/^Sbjct:\s+\d+\s+[A-Z,\-,\*]+\s+(\d+)$/i) { |
| 285 |
|
|
$h3=$1; |
| 286 |
|
|
} |
| 287 |
|
|
else { |
| 288 |
|
|
&err_parse("Cannot find end of Sbjct: line for alignment (query=$query, subj=$hit)"); |
| 289 |
|
|
} |
| 290 |
|
|
$fend = $bytes_read; |
| 291 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 292 |
|
|
} |
| 293 |
|
|
$fend=$bytes_read; |
| 294 |
|
|
if ($curgrp) { |
| 295 |
|
|
push(@grpdata, [$q5, $q3, $h5, $h3, $pid, $psim, $score, $frame, $ids, $sims, $bcount]); |
| 296 |
|
|
#$grphspcount++; |
| 297 |
|
|
#$grpids+=$ids; |
| 298 |
|
|
#$grpsims+=$sims; |
| 299 |
|
|
#$grplen+=$bcount; |
| 300 |
|
|
} |
| 301 |
|
|
else { |
| 302 |
|
|
&$writeline; |
| 303 |
|
|
} |
| 304 |
|
|
next; #skip to next HSP, as we might have read it |
| 305 |
|
|
} # alignment parsing |
| 306 |
|
|
} |
| 307 |
|
|
SKIP_EMPTYLN: |
| 308 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 309 |
|
|
&err_parse("--Timeout trying to find all HSPs for query $query and subj $hit") |
| 310 |
|
|
if ($.>$savepos+5000); |
| 311 |
|
|
} |
| 312 |
|
|
&err_parse("--Cannot get hit data for $hit") |
| 313 |
|
|
if ($score==0 || $h3==0); |
| 314 |
|
|
|
| 315 |
|
|
} #--Subject entry |
| 316 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 317 |
|
|
} |
| 318 |
|
|
NEXT_QUERY: |
| 319 |
|
|
do { |
| 320 |
|
|
$_=<STDIN>;$bytes_read+=length($_); |
| 321 |
|
|
last if eof(STDIN); |
| 322 |
|
|
#} until /^\x0C/; # record separator for multi-fasta wu-blast! |
| 323 |
|
|
} until (/^Query=/); |
| 324 |
|
|
goto NEXT_QRYBLOCK; |
| 325 |
|
|
} |
| 326 |
|
|
|
| 327 |
|
|
exit(0); |
| 328 |
|
|
|
| 329 |
|
|
#========================================== |
| 330 |
|
|
#to speed it up, use the variables directly from main program |
| 331 |
|
|
#(not very elegant, but more efficient) |
| 332 |
|
|
|
| 333 |
|
|
sub err_parse { |
| 334 |
|
|
die "Parsing Error at line $.:\n$_\n$_[0]\n"; |
| 335 |
|
|
} |
| 336 |
|
|
|
| 337 |
|
|
sub write_noflt { |
| 338 |
|
|
$hdefline =~ s/^>/\x01/; |
| 339 |
|
|
$hdefline =~ s/ >([\w\|]{6,})/\x01$1/g; |
| 340 |
|
|
my $xtra=$_[0]?"\t$_[0]":""; |
| 341 |
|
|
print join("\t", |
| 342 |
|
|
# 0 1 2 3 4 5 6 7 8 9 10 11 |
| 343 |
|
|
($query, $date, $qlen, $method, $searchdb, $hit, $q5, $q3, $h5, $h3, $pid, $psim, |
| 344 |
|
|
# 12 13 14 15 16 17 18 19 |
| 345 |
|
|
$score, $fpos, $fend, $hdefline, $frame, $strand, $hitlen, $evalue))."$xtra\n"; |
| 346 |
|
|
}; #clean write |
| 347 |
|
|
|
| 348 |
|
|
sub write_flt { |
| 349 |
|
|
#my ($minhlen, $maxoh, $minpid, $minsc, $maxE) |
| 350 |
|
|
my ($minmatch, $maxovhang)=($is_lenperc?(($minhlen*$qlen)/100):$minhlen, |
| 351 |
|
|
$is_ovperc?(($maxoh*$qlen)/100):$maxoh); |
| 352 |
|
|
|
| 353 |
|
|
my ($q_match, $q_beg_over, $q_end_over)=($q5 < $q3)? |
| 354 |
|
|
(($q3 - $q5) + 1, $q5 - 1, $qlen - $q3) : |
| 355 |
|
|
(($q5 - $q3) + 1, $qlen - $q5, $q3 - 1); |
| 356 |
|
|
my ($h_match, $h_beg_over, $h_end_over)= ($h5 < $h3)? |
| 357 |
|
|
( ($h3 - $h5) + 1, $h5 - 1, $hitlen - $h3) : |
| 358 |
|
|
( ($h5 - $h3) + 1, $hitlen - $h5, $h3 - 1); |
| 359 |
|
|
&write_noflt($_[0]) |
| 360 |
|
|
if (($query ne $hit) && |
| 361 |
|
|
(($h_match >=$minmatch) || ($q_match >= $minmatch)) |
| 362 |
|
|
&& (($q_beg_over <= $maxovhang) || ($h_beg_over <= $maxovhang)) |
| 363 |
|
|
&& (($q_end_over <= $maxovhang) || ($h_end_over <= $maxovhang)) |
| 364 |
|
|
&& ($pid>=$minpid) |
| 365 |
|
|
&& ($psim>=$minpsim) |
| 366 |
|
|
&& ($score>=$minscore) |
| 367 |
|
|
&& ($evalue<=$maxE)) |
| 368 |
|
|
} |
| 369 |
|
|
|
| 370 |
|
|
sub write_group { |
| 371 |
|
|
die "Error: write_group called with empty \@grpdata\n" |
| 372 |
|
|
if (@grpdata==0); |
| 373 |
|
|
my (@q5, @q3, @h5, @h3); |
| 374 |
|
|
@grpdata = sort { $main::a->[0] <=> $main::b->[0] } @grpdata; |
| 375 |
|
|
my @details; |
| 376 |
|
|
$score=0; |
| 377 |
|
|
my $hitcov; #non redundant hit coverage -- |
| 378 |
|
|
my @hitsegs; # non-overlapping hit segments |
| 379 |
|
|
my $frm; |
| 380 |
|
|
my ($grplen, $grpids, $grpsims)=(0,0,0); |
| 381 |
|
|
# also recompute $psim, $pid, $q5, $q3, $h5, $h3 |
| 382 |
|
|
# because psim/pid filter may have discarded some HSPs |
| 383 |
|
|
#my ($sumpid, $sumpsim, $sumlen, $minq, $maxq, $minh, $maxh); |
| 384 |
|
|
foreach my $d (@grpdata) { |
| 385 |
|
|
#$d=[$q5, $q3, $h5, $h3, $pid, $psim, $score, $frame, $numids, $numsims, $numbases] |
| 386 |
|
|
# 0 1 2 3 4 5 6 7 8 9 10 |
| 387 |
|
|
#print STDERR "seg $$d[0]-$$d[1] : $$d[2]-$$d[3]\n"; |
| 388 |
|
|
next unless ($$d[4]>=$minpid && $$d[5]>=$minpsim); |
| 389 |
|
|
$grpids += $$d[8]; |
| 390 |
|
|
$grpsims += $$d[9]; |
| 391 |
|
|
$grplen += $$d[10]; |
| 392 |
|
|
my ($seg_start, $seg_end)=($$d[1]<$$d[0])?($$d[1],$$d[0]) :($$d[0],$$d[1]); |
| 393 |
|
|
my $seglen=$seg_end-$seg_start+1; |
| 394 |
|
|
#$sumpid+=$seglen*$$d[4]; |
| 395 |
|
|
#$sumpsim+=$seglen*$$d[5]; |
| 396 |
|
|
#$sumlen+=$seglen; |
| 397 |
|
|
#$minq=$seg_start if ($minq==0 || $minq>$seg_start); |
| 398 |
|
|
#$maxq=$seg_end if ($maxq==0 || $maxq<$seg_end); |
| 399 |
|
|
|
| 400 |
|
|
push(@q5,$$d[0]);push(@q3,$$d[1]); |
| 401 |
|
|
push(@h5,$$d[2]);push(@h3,$$d[3]); |
| 402 |
|
|
if ($$d[2]>$$d[3]) { |
| 403 |
|
|
#$minh=$$d[3] if ($minh==0 || $minh>$$d[3]); |
| 404 |
|
|
#$maxh=$$d[2] if ($maxh==0 || $maxh<$$d[2]); |
| 405 |
|
|
push(@hitsegs,[$$d[3],$$d[2]] ) |
| 406 |
|
|
} |
| 407 |
|
|
else { |
| 408 |
|
|
#$minh=$$d[3] if ($minh==0 || $minh>$$d[3]); |
| 409 |
|
|
#$maxh=$$d[2] if ($maxh==0 || $maxh<$$d[2]); |
| 410 |
|
|
push(@hitsegs,[$$d[2],$$d[3]]); |
| 411 |
|
|
} |
| 412 |
|
|
push(@details, "$$d[0]-$$d[1]\:$$d[2]-$$d[3]|$$d[5]"); |
| 413 |
|
|
$score+=$$d[6]; |
| 414 |
|
|
$frm=$$d[7] unless $frm; |
| 415 |
|
|
$frame='*' if ($$d[7] ne $frm); |
| 416 |
|
|
} #for each HSP |
| 417 |
|
|
#if ($sumlen>0) { |
| 418 |
|
|
# $psim=$sumpsim/$sumlen; |
| 419 |
|
|
# $pid=$sumpid/$sumlen; |
| 420 |
|
|
# } |
| 421 |
|
|
goto GO_OUT if ($grplen==0); |
| 422 |
|
|
$pid=sprintf('%d',$grpids*100/$grplen); |
| 423 |
|
|
$psim=sprintf('%d',$grpsims*100/$grplen); |
| 424 |
|
|
|
| 425 |
|
|
goto GO_OUT if ($psim<$minpsim || $pid <$minpid); |
| 426 |
|
|
|
| 427 |
|
|
$hitcov=&mergeSegs(\@hitsegs); |
| 428 |
|
|
|
| 429 |
|
|
#@q5 = sort { $main::a <=> $main::b } @q5; |
| 430 |
|
|
#@q3 = sort { $main::a <=> $main::b } @q3; |
| 431 |
|
|
#@h5 = sort { $main::a <=> $main::b } @h5; |
| 432 |
|
|
#@h3 = sort { $main::a <=> $main::b } @h3; |
| 433 |
|
|
#determine actual (non-overlapping) coverage of the hit |
| 434 |
|
|
($q5, $q3, $h5, $h3)=($q5[0], $q3[-1], $h5[0], $h3[-1]); |
| 435 |
|
|
if ($strand eq 'Minus' && ($q5[0]>$q3[0])) { # q5>q3 unless protein search |
| 436 |
|
|
($q5, $q3, $h5, $h3)=($q5[-1], $q3[0], $h5[0], $h3[-1]); |
| 437 |
|
|
@details=reverse(@details); |
| 438 |
|
|
} |
| 439 |
|
|
|
| 440 |
|
|
&$writeline("$hitcov\t".join('~',@details)); |
| 441 |
|
|
#-- clear |
| 442 |
|
|
GO_OUT: |
| 443 |
|
|
@grpdata=(); |
| 444 |
|
|
undef($curgrp); |
| 445 |
|
|
#$grplen=0;$grpids=0;$grpsims=0;$grphspcount=0; |
| 446 |
|
|
} |
| 447 |
|
|
|
| 448 |
|
|
sub mergeSegs { |
| 449 |
|
|
my $a=shift; # ref to list of [h5, h3] |
| 450 |
|
|
my $hcov=0; |
| 451 |
|
|
OVLCHECK: |
| 452 |
|
|
for (my $i=0;$i<(@$a-1); $i++) { |
| 453 |
|
|
for (my $j=$i+1;$j<@$a; $j++) { |
| 454 |
|
|
my ($li, $ri, $lj, $rj)= |
| 455 |
|
|
($$a[$i]->[0], $$a[$i]->[1], $$a[$j]->[0], $$a[$j]->[1]); |
| 456 |
|
|
if ($lj<=$ri && $lj>=$li) { |
| 457 |
|
|
$$a[$i]->[0]=$li; |
| 458 |
|
|
$$a[$i]->[1]=($rj>$ri) ? $rj : $ri; |
| 459 |
|
|
splice(@$a, $j,1); |
| 460 |
|
|
goto OVLCHECK; |
| 461 |
|
|
} |
| 462 |
|
|
elsif ($li<=$rj && $li>=$lj) { |
| 463 |
|
|
$$a[$i]->[0]=$lj; |
| 464 |
|
|
$$a[$i]->[1]=($rj>$ri) ? $rj : $ri; |
| 465 |
|
|
splice(@$a, $j,1); |
| 466 |
|
|
goto OVLCHECK; |
| 467 |
|
|
} |
| 468 |
|
|
} # for j |
| 469 |
|
|
} # for i |
| 470 |
|
|
foreach my $d (@$a) { |
| 471 |
|
|
$hcov+=($$d[1]-$$d[0]+1); |
| 472 |
|
|
} |
| 473 |
|
|
return $hcov; |
| 474 |
|
|
} |
