| 1 |
#!/usr/bin/perl |
| 2 |
use strict; |
| 3 |
use Getopt::Std; |
| 4 |
use FindBin;use lib $FindBin::Bin; |
| 5 |
|
| 6 |
my $usage = q/Usage: |
| 7 |
cuffdiff_loc_check.pl [-c <min_cov>] [-f <min_fpkm>] [-l <labels>] [-b <bamfiles>] \ |
| 8 |
<gene_exp.diff> <cuffcompare.loci> <cuffcompare.tracking> <transfrags1.ifa> <transfrags2.ifa> |
| 9 |
|
| 10 |
(assume that transfrags?.ifa files were created using gff2iit -XTA ) |
| 11 |
If BAM files are provided for specificity validation, they must be indexed. |
| 12 |
/; |
| 13 |
umask 0002; |
| 14 |
getopts('b:f:c:l:o:') || die($usage."\n"); |
| 15 |
my $outfile=$Getopt::Std::opt_o; |
| 16 |
my $lnames=$Getopt::Std::opt_l; |
| 17 |
my $bamfiles=$Getopt::Std::opt_b; |
| 18 |
my $MIN_COV=$Getopt::Std::opt_c || 2; |
| 19 |
my $MIN_FPKM=$Getopt::Std::opt_f || 2; |
| 20 |
if ($outfile) { |
| 21 |
open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n"); |
| 22 |
select(OUTF); |
| 23 |
} |
| 24 |
|
| 25 |
my @labels=split(/[\,\s\;]/,$lnames) if $lnames; |
| 26 |
my @bam=split(/[\,\s\;]/,$bamfiles) if $bamfiles; |
| 27 |
foreach my $bfile (@bam) { |
| 28 |
die("Error: file $bfile not found!\n") unless -f $bfile; |
| 29 |
die("Error: bam index $bfile.bai not found!\n") unless -f $bfile.'.bai'; |
| 30 |
} |
| 31 |
unless ($labels[1]) { |
| 32 |
@labels=('Exp1','Exp2'); |
| 33 |
} |
| 34 |
# -- |
| 35 |
my ($gfexp, $lfile, $kfile, $tfile1, $tfile2)=@ARGV; |
| 36 |
die($usage) unless $tfile2; |
| 37 |
# 0 1 2 3 4 5 6 |
| 38 |
my %locdiff; # XLOC_id => [chr, c_start, c_end, fpkm_1, fpkm_2, p-value, is_significant] |
| 39 |
# cuffdiff loci data |
| 40 |
|
| 41 |
# 0 1 2 3 4 5 6 |
| 42 |
my %loci; # locus_id => [chr, strand, start, end, [refs], [@ts1], [@ts2]] |
| 43 |
my %t1; # t_id => [$loc, ref, code, FPKM, cov, covlen, $exd] |
| 44 |
my %t2; # t_id => [$loc, ref, code, FPKM, cov, covlen, $exd] |
| 45 |
# 0 1 2 3 4 5 6 |
| 46 |
open(DFILE, $gfexp) || die("Error opening file $gfexp\n"); |
| 47 |
|
| 48 |
while (<DFILE>) { |
| 49 |
next unless m/^XLOC/; |
| 50 |
chomp; |
| 51 |
my @t=split(/\t/); |
| 52 |
#next unless $t[-1] eq 'yes'; |
| 53 |
my $is_significant = ($t[-1] eq 'yes'); |
| 54 |
my ($chr, $cstart, $cend)=($t[2]=~m/^([^\:]+):(\d+)\-(\d+)/); |
| 55 |
($cstart, $cend)=($cend,$cstart) if $cend<$cstart; |
| 56 |
die("Error: no chromosome coordinates parsed ($_)!\n") unless $cstart && $cend; |
| 57 |
$locdiff{$t[0]}=[$chr, $cstart, $cend, $t[6], $t[7], $t[10], $is_significant]; |
| 58 |
} |
| 59 |
close(DFILE); |
| 60 |
die("Error: no gene diff data loaded from $lfile\n") unless keys(%locdiff)>0; |
| 61 |
print STDERR "Loading $gfexp file..\n"; |
| 62 |
|
| 63 |
open(LFILE, $lfile) || die("Error opening file $lfile!\n"); |
| 64 |
print STDERR "Loading loci file..\n"; |
| 65 |
while (<LFILE>) { |
| 66 |
chomp; |
| 67 |
my @t=split(/\t/); |
| 68 |
my ($chr, $strand, $cstart, $cend)= |
| 69 |
($t[1]=~m/^([^\[]+)\[([\+\-\.])\](\d+)\-(\d+)/); |
| 70 |
my @refs=split(/\,/,$t[2]) if $t[2] ne '-'; |
| 71 |
my @ts1=split(/\,/,$t[3]) if $t[3] ne '-'; |
| 72 |
my @ts2=split(/\,/,$t[4]) if $t[4] ne '-'; |
| 73 |
$loci{$t[0]}=[$chr, $strand, $cstart, $cend, [@refs], [@ts1], [@ts2]]; |
| 74 |
} |
| 75 |
close(LFILE); |
| 76 |
print STDERR "Loading tracking file..\n"; |
| 77 |
open(KFILE, $kfile) || die("Error opening file $kfile!\n"); |
| 78 |
while (<KFILE>) { |
| 79 |
chomp; |
| 80 |
my @t=split(/\t/); |
| 81 |
#next unless ($t[2] eq '-' && $t[3] ne 'r') || $t[3]=~tr/.p\-ucj=eoi//; |
| 82 |
next if $t[3]=~tr/xs//; # tr/xsr// |
| 83 |
$t[2]='' if $t[2]=~m/^\s*\-\s*$/; # reference transcript |
| 84 |
unless ($t[4] eq '-') { |
| 85 |
my @q=split(/\|/,$t[4]); |
| 86 |
$t1{$q[1]}=[$t[1], $t[2], $t[3], $q[3], $q[6], 0, []]; |
| 87 |
} |
| 88 |
unless ($t[5] eq '-') { |
| 89 |
my @q=split(/\|/,$t[5]); |
| 90 |
$t2{$q[1]}=[$t[1], $t[2], $t[3], $q[3], $q[6], 0, []]; |
| 91 |
} |
| 92 |
if ($t[2]) { |
| 93 |
# update ref list in loci as sometimes it misses some o: and i: genes (!) |
| 94 |
my $ld=$loci{$t[1]}; |
| 95 |
die("Error: locus $t[1] not found in the loci file but found in tracking file!\n") unless $ld; |
| 96 |
my %g; |
| 97 |
if (@{$ld->[4]}>0) { |
| 98 |
@g{@{$ld->[4]}}=(); |
| 99 |
} |
| 100 |
$g{$t[2]}=1; |
| 101 |
$ld->[4]=[(keys(%g))]; |
| 102 |
} |
| 103 |
} |
| 104 |
close(KFILE); |
| 105 |
|
| 106 |
print STDERR "Loading transcript file 1 ..\n"; |
| 107 |
open(TFILE, $tfile1) || die("Error opening file $tfile1!\n"); |
| 108 |
my $t; #current transcript |
| 109 |
while (<TFILE>) { |
| 110 |
chomp; |
| 111 |
if (m/^>(\S+)\s+([^:]+):(\d+)\.\.(\d+)/) { |
| 112 |
my ($tid, $chr, $cstart, $cend)=($1,$2,$3,$4); |
| 113 |
my $td=$t1{$tid}; |
| 114 |
next unless $td; |
| 115 |
my ($covlen, $numexons, @ex); |
| 116 |
while (<TFILE>) { |
| 117 |
chomp; |
| 118 |
if (m/^\d+\-/) { |
| 119 |
@ex=map { [split(/\-/)] } (split(/\,/)); |
| 120 |
$numexons=scalar(@ex); |
| 121 |
map { $covlen+=$_->[1]-$_->[0]+1 } @ex; |
| 122 |
next; |
| 123 |
} |
| 124 |
if (m/^i:/) { |
| 125 |
my ($fpkm)=(m/FPKM=([^;]+)/); |
| 126 |
my ($cov)=(m/cov=([^;]+)/); |
| 127 |
($$td[3], $$td[4], $$td[5], $$td[6])= |
| 128 |
($fpkm, $cov, $covlen, [@ex]); |
| 129 |
# -- debug: |
| 130 |
#print STDERR "..saving ".scalar(@ex)." exons \n"; |
| 131 |
#my ($loc, $ref, $code, $fpk, $cc, $covl, $x)=@$td; |
| 132 |
#my $numex=scalar(@$x); |
| 133 |
#die "$tid: $loc $ref\[$code\] cov=$cc covlen=$covl numexons=$numex\n"; |
| 134 |
#-- |
| 135 |
last; |
| 136 |
} |
| 137 |
} #while record body |
| 138 |
} #header found |
| 139 |
} |
| 140 |
close(TFILE); |
| 141 |
print STDERR "Loading transcript file 2 ..\n"; |
| 142 |
open(TFILE, $tfile2) || die("Error opening file $tfile2!\n"); |
| 143 |
while (<TFILE>) { |
| 144 |
if (m/^>(\S+)\s+([^:]+):(\d+)\.\.(\d+)/) { |
| 145 |
my ($tid, $chr, $cstart, $cend)=($1,$2,$3,$4); |
| 146 |
my $td=$t2{$tid}; |
| 147 |
next unless $td; |
| 148 |
my ($covlen, $numexons, @ex); |
| 149 |
while (<TFILE>) { |
| 150 |
chomp; |
| 151 |
if (m/^\d+\-/) { |
| 152 |
@ex=map { [split(/\-/)] } (split(/\,/)); |
| 153 |
$numexons=scalar(@ex); |
| 154 |
map { $covlen+=$_->[1]-$_[0]+1 } @ex; |
| 155 |
next; |
| 156 |
} |
| 157 |
if (m/^i:/) { |
| 158 |
my ($fpkm)=(m/FPKM=([^;]+)/); |
| 159 |
my ($cov)=(m/cov=([^;]+)/); |
| 160 |
($$td[3], $$td[4], $$td[5], $$td[6])= |
| 161 |
($fpkm, $cov, $covlen, [@ex]); |
| 162 |
last; |
| 163 |
} |
| 164 |
} #while record body |
| 165 |
} #header found |
| 166 |
} |
| 167 |
close(TFILE); |
| 168 |
print STDERR "Now compute coverage and print results.\n"; |
| 169 |
my $maxcoverage=0; |
| 170 |
my $max_fpkm=0; |
| 171 |
my ($maxcov_loc, $maxfpkm_loc); |
| 172 |
# now finally check each locus |
| 173 |
while (my ($l,$ld)=each(%loci)) { |
| 174 |
my ($chr, $strand, $cstart, $cend, $rrefs, $rt1, $rt2)=@$ld; |
| 175 |
my @refs=@$rrefs; |
| 176 |
# check if this "prediction" is confirmed by cuffdiff |
| 177 |
my $tlist1 = (@$rt1>0) ? join(',',@$rt1) : '-'; |
| 178 |
my $tlist2 = (@$rt2>0) ? join(',',@$rt2) : '-'; |
| 179 |
#merge all exons of @ts1 |
| 180 |
#my $debug=("$chr:$cstart-$cend" eq 'chr3:41288638-41289701'); |
| 181 |
my ($mcov1, $mcovlen1, $maxcovlen1, $maxexons1, $maxfpkm1, $dcode1)=getAvgCov($rt1, \%t1); |
| 182 |
my ($mcov2, $mcovlen2, $maxcovlen2, $maxexons2, $maxfpkm2, $dcode2)=getAvgCov($rt2, \%t2); |
| 183 |
($maxcoverage, $maxcov_loc)=($mcov1,$l) if $mcov1>$maxcoverage; |
| 184 |
($maxcoverage, $maxcov_loc)=($mcov2,$l) if $mcov2>$maxcoverage; |
| 185 |
($max_fpkm, $maxfpkm_loc)=($maxfpkm1,$l) if $maxfpkm1>$max_fpkm; |
| 186 |
($max_fpkm, $maxfpkm_loc)=($maxfpkm2,$l) if $maxfpkm2>$max_fpkm; |
| 187 |
next if ($maxcovlen1==0 && $maxcovlen2==0); |
| 188 |
my $dcode=$dcode1.$dcode2; |
| 189 |
my $refcode; |
| 190 |
if (length($dcode)==2 && $dcode1 eq $dcode2 && $dcode1 ne '.') { |
| 191 |
$refcode=$dcode1.':'; # prefix for reference column |
| 192 |
} |
| 193 |
elsif (length($dcode)==1 && $dcode ne '.') { |
| 194 |
$refcode=$dcode.':'; |
| 195 |
} |
| 196 |
#print STDERR "$l has uniform refcode: $refcode\n" if $refcode; |
| 197 |
my $reflist = (@refs>0) ? join(',',@refs) : '-'; |
| 198 |
#-- discard single exon "novel" genes for now |
| 199 |
# (though they could be some relevant non-coding RNAs..) |
| 200 |
next if ($reflist eq '-' || $strand eq '.') && $maxexons1<2 && $maxexons2<2; |
| 201 |
my $status; |
| 202 |
my $pnovel=($reflist eq '-' || $refcode eq 'i:' || $refcode eq 'p:'); # has potential to be a novel gene |
| 203 |
my $isnovel=0; |
| 204 |
#my ($fpkm1, $fpkm2, $pval)=('-','-','-'); |
| 205 |
#if ($pnovel && ($mcov1+$mcov2>=$MIN_COV)) { |
| 206 |
# $status='novel'; |
| 207 |
# $isnovel=1; |
| 208 |
# } |
| 209 |
|
| 210 |
# if not "novel", discard it unless |
| 211 |
# cuffdiff found a significant change in expression level |
| 212 |
# in a region that is consistent with this locus |
| 213 |
my $diffd=$locdiff{$l}; |
| 214 |
my ($dchr, $dstart, $dend, $fpkm1, $fpkm2, $pval, $is_significant); |
| 215 |
unless ($diffd) { |
| 216 |
#print STDERR "Warning: $l not found in $gfexp file, skipping..\n"; |
| 217 |
#next; |
| 218 |
($fpkm1, $fpkm2)=($maxfpkm1, $maxfpkm2); # may not be reliable.. |
| 219 |
next if ($fpkm1+$fpkm2<$MIN_FPKM && $mcov1+$mcov2<$MIN_COV); |
| 220 |
if ($pnovel) { |
| 221 |
$status='novel'; |
| 222 |
$isnovel=1; |
| 223 |
} |
| 224 |
} |
| 225 |
else { |
| 226 |
#next unless ($isnovel || $diffd); |
| 227 |
#if ($diffd) { |
| 228 |
($dchr, $dstart, $dend, $fpkm1, $fpkm2, $pval, $is_significant)=@$diffd; |
| 229 |
#($fpkm1, $fpkm2, $pval)=($v1,$v2, $dp); |
| 230 |
if ($chr ne $dchr) { |
| 231 |
print STDERR "Warning: $l not found on the same chromosome in $gfexp file?\n"; |
| 232 |
next; |
| 233 |
} |
| 234 |
if ($pnovel && ($is_significant || ($fpkm1+$fpkm2>=$MIN_FPKM && $mcov1+$mcov2>=$MIN_COV))) { |
| 235 |
$status='novel'; |
| 236 |
$isnovel=1; |
| 237 |
} |
| 238 |
next unless ($isnovel || $is_significant); |
| 239 |
my $loclen=($cend-$cstart+1); |
| 240 |
my $dloclen=($dend-$dstart+1); |
| 241 |
# check if the quantification interval can be trusted |
| 242 |
if (!$isnovel && ($dloclen-$loclen>400 || $dend-$cend>200 || $cstart-$dstart>200)) { |
| 243 |
# bad quantification interval |
| 244 |
$is_significant=0; |
| 245 |
#next if ($dloclen-$loclen>400 || $dend-$cend>200 || $cstart-$dstart>200); |
| 246 |
} |
| 247 |
} |
| 248 |
my $libtag=$isnovel ? '-novel' : '-specific'; |
| 249 |
if ($tlist1 eq '-' && $fpkm1==0 && ($is_significant || ($fpkm2>=$MIN_FPKM && $mcov2>=$MIN_COV))) { |
| 250 |
if ($bamfiles && !bamHasOvl($chr, $cstart, $cend, $rt2, \%t2, $bam[0])) { |
| 251 |
$status=$labels[1].$libtag; |
| 252 |
} |
| 253 |
} |
| 254 |
elsif ($tlist2 eq '-' && $fpkm2==0 && ($is_significant || ($fpkm1>=$MIN_FPKM && $mcov1>=$MIN_COV))) { |
| 255 |
if ($bamfiles && !bamHasOvl($chr, $cstart, $cend, $rt1, \%t1, $bam[1])) { |
| 256 |
$status=$labels[0].$libtag; |
| 257 |
} |
| 258 |
} |
| 259 |
unless ($status) { |
| 260 |
next if $refcode eq 'o:'; #discard dubious overlaps |
| 261 |
if ($mcov1+$mcov2>0 && $fpkm1+$fpkm2>0) { |
| 262 |
if ($is_significant) { |
| 263 |
$status = ($fpkm1>$fpkm2) ? $labels[0].'-inc_cdiff' : $labels[1].'-inc_cdiff'; |
| 264 |
} |
| 265 |
#else { |
| 266 |
my $covratio=signedRatio($mcov1, $mcov2); #decrease ratio |
| 267 |
my $covlenratio=signedRatio($mcovlen1,$mcovlen2); #decrease ratio, negative if there is an increase |
| 268 |
if (($covlenratio>=2 && $covratio>=3) || (abs($covlenratio)<2 && $covratio>4)) { |
| 269 |
updStat(\$status, $labels[0],'inc_cov') if $fpkm1>$fpkm2; |
| 270 |
} |
| 271 |
elsif (($covlenratio<=-2 && $covratio<=-3) || (abs($covlenratio)<2 && $covratio<-4)) { |
| 272 |
updStat(\$status, $labels[1],'inc_cov') if $fpkm2>$fpkm1; |
| 273 |
} |
| 274 |
elsif (abs($covlenratio)<2) { |
| 275 |
if ($fpkm1>2*$fpkm2 && $mcov1>$mcov2*2) { updStat(\$status, $labels[0],'inc_cov') } |
| 276 |
elsif ($fpkm2>2*$fpkm1 && $mcov2>$mcov1*2) { updStat(\$status, $labels[1],'inc_cov')} |
| 277 |
} |
| 278 |
# } |
| 279 |
} |
| 280 |
} |
| 281 |
next unless $status; |
| 282 |
|
| 283 |
print join("\t",$chr, $strand, $cstart, $cend, "$chr\:$cstart-$cend", $l, $status, |
| 284 |
$refcode.$reflist, join('|', $mcov1, $maxexons1, $maxfpkm1, $tlist1), |
| 285 |
join('|', $mcov2, $maxexons2, $maxfpkm2, $tlist2), $fpkm1, $fpkm2, $pval)."\n"; |
| 286 |
} # for each locus |
| 287 |
print STDERR "Max. coverage found: $maxcoverage for locus $maxcov_loc\n"; |
| 288 |
print STDERR "Max. FPKM found: $max_fpkm for locus $maxfpkm_loc\n"; |
| 289 |
# -- |
| 290 |
if ($outfile) { |
| 291 |
select(STDOUT); |
| 292 |
close(OUTF); |
| 293 |
} |
| 294 |
|
| 295 |
#************ Subroutines ************** |
| 296 |
sub updStat { |
| 297 |
my ($rstatus, $lib, $tag)=@_; |
| 298 |
if (length($$rstatus)==0) { |
| 299 |
$$rstatus=$lib.'-'.$tag; |
| 300 |
} |
| 301 |
else { |
| 302 |
$$rstatus.='-'.$tag; |
| 303 |
} |
| 304 |
} |
| 305 |
|
| 306 |
sub getAvgCov { |
| 307 |
my ($tlst, $th)=@_; |
| 308 |
my $maxcovlen=0; |
| 309 |
my $tcov=0; |
| 310 |
my $tcovlen=0; |
| 311 |
my $tcount=0; |
| 312 |
my $maxexons=0; |
| 313 |
my %codes; |
| 314 |
my $domcode=''; |
| 315 |
my $maxcov=0; |
| 316 |
my $maxfpkm=0; |
| 317 |
foreach my $t (@$tlst) { |
| 318 |
my $thd=$th->{$t}; |
| 319 |
next unless ($thd); |
| 320 |
my ($loc, $ref, $code, $fpkm, $cov, $covlen, $ex)=@$thd; |
| 321 |
$domcode='.' if $ref; |
| 322 |
$codes{$code}++ unless $code=~tr/u.\-r//; |
| 323 |
my $numexons=scalar(@$ex); |
| 324 |
#die "$t: $loc $ref\[$code\] covlen=$covlen, exons=$numexons\n"; |
| 325 |
|
| 326 |
$maxcovlen=$covlen if ($maxcovlen<$covlen); |
| 327 |
$maxexons=$numexons if $maxexons<$numexons; |
| 328 |
$tcovlen+=$covlen; |
| 329 |
$tcount++; |
| 330 |
$maxcov=$cov if ($cov>$maxcov); |
| 331 |
$maxfpkm=$fpkm if $fpkm>$maxfpkm; |
| 332 |
$tcov+=$cov; |
| 333 |
} |
| 334 |
my @classes=keys(%codes); |
| 335 |
$domcode=$classes[0] if (@classes==1); |
| 336 |
#my ($mcov, $mcovlen)=($tcount==0) ? ( 0, 0 ) : |
| 337 |
# (sprintf('%.3f', $tcov/$tcount), sprintf('%d', $tcovlen/$tcount)); |
| 338 |
my ($mcov, $mcovlen)=($tcount==0) ? ( 0, 0 ) : (sprintf('%.2f',$maxcov), $maxcovlen); |
| 339 |
return ($mcov, $mcovlen, $maxcovlen, $maxexons, sprintf('%.2f',$maxfpkm), $domcode); |
| 340 |
} |
| 341 |
|
| 342 |
sub signedRatio { |
| 343 |
# decrease ratio from $a to $b (negative if there is an increase) |
| 344 |
my ($a,$b)=@_; |
| 345 |
if ($a>$b) { return ($b==0)?1000 : $a/$b } |
| 346 |
elsif ($a==$b) { |
| 347 |
return ($a==0) ? 0 : 1 |
| 348 |
} |
| 349 |
else { return ($a==0) ? -1000 : -$b/$a } |
| 350 |
} |
| 351 |
|
| 352 |
sub bamHasOvl { |
| 353 |
my ($chr, $cstart, $cend, $tlst, $th, $fbam)=@_; |
| 354 |
|
| 355 |
my $range=$chr.':'.$cstart.'-'.$cend; |
| 356 |
#my $debug=($range eq 'chr3:41288638-41289701'); |
| 357 |
my @exons; |
| 358 |
foreach my $t (@$tlst) { |
| 359 |
my $thd=$th->{$t}; |
| 360 |
next unless $thd; |
| 361 |
my ($loc, $ref, $code, $fpkm, $cov, $covlen, $ex)=@$thd; |
| 362 |
push(@exons, @$ex); |
| 363 |
} #gather exons from all transcripts in this locus |
| 364 |
return 0 unless @exons>0; |
| 365 |
my @sorted_exons= sort { $a->[0]<=>$b->[0] } @exons; |
| 366 |
open(SAMPIPE, "samtools view $fbam $range |") || die ("Error opening samtools pipe ($!)\n"); |
| 367 |
while(<SAMPIPE>) { |
| 368 |
my $samline=$_; |
| 369 |
chomp; |
| 370 |
my ($qname, $flags, $gseq, $pos, $mapq, $cigar, $rnext, $pnext, $tlen, $seq, $qual, @extra)= |
| 371 |
split(/\t/); |
| 372 |
# if ($strand) { |
| 373 |
# my $mstrand= (($flags & 0x10)==0) ? '+' : '-'; |
| 374 |
# next if $mstrand ne $strand; |
| 375 |
# } |
| 376 |
#now extract the CIGAR segments |
| 377 |
my @cigdata=($cigar=~m/(\d+[A-Z,=])/g); |
| 378 |
my ($mstart, $mend); |
| 379 |
my $hasOvl=0; |
| 380 |
my $curpos=$pos; |
| 381 |
$mstart=$pos; |
| 382 |
foreach my $cd (@cigdata) { |
| 383 |
my $code=chop($cd); |
| 384 |
if ($code eq 'N') { #gap |
| 385 |
#process previous interval |
| 386 |
if ($mend && checkOverlap($mstart, $mend, \@sorted_exons)) { |
| 387 |
$hasOvl=1; |
| 388 |
last; |
| 389 |
} |
| 390 |
$mstart=$curpos+$cd; |
| 391 |
$mend=undef; |
| 392 |
next; |
| 393 |
} |
| 394 |
$mend=$curpos+$cd-1; |
| 395 |
$curpos+=$cd; |
| 396 |
} |
| 397 |
unless ($hasOvl) { #check the last interval |
| 398 |
if ($mend && checkOverlap($mstart, $mend, \@sorted_exons)) { |
| 399 |
$hasOvl=1; |
| 400 |
} |
| 401 |
} |
| 402 |
if ($hasOvl) { |
| 403 |
close(SAMPIPE); |
| 404 |
return 1; |
| 405 |
} |
| 406 |
} # while <SAMPIPE> |
| 407 |
close(SAMPIPE); |
| 408 |
return 0; |
| 409 |
} |
| 410 |
|
| 411 |
|
| 412 |
|
| 413 |
sub checkOverlap { |
| 414 |
my ($a, $b, $rx)=@_; |
| 415 |
return 0 if ($a>$$rx[-1]->[1] || $b<$$rx[0]->[0]); # not overlapping the whole exon chain |
| 416 |
foreach my $x (@$rx) { |
| 417 |
#return (start<=d->end && end>=d->start) |
| 418 |
return 1 if ($a<=$$x[1] && $b>=$$x[0]); |
| 419 |
return 0 if $b<$$x[0]; #@$rx is sorted |
| 420 |
} |
| 421 |
} |
