| 1 |
#!/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 |
} |
