File: changepointCaller.pm

package info (click to toggle)
alien-hunter 1.7-8
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, bullseye, sid
  • size: 240 kB
  • sloc: perl: 962; java: 358; sh: 32; makefile: 9
file content (126 lines) | stat: -rw-r--r-- 2,647 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
=head1 NAME

changepointCaller

=head1 SYNOPSIS

changepointCaller - it calls the PAI_scripts::changePoint

=head1 AUTHOR

George Vernikos <gsv(at)sanger.ac.uk>


=head1 LICENSE

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

=cut

use FindBin;
use lib $FindBin::Bin;

package changepointCaller;
use PAI_scripts::changePoint;
use Exporter;
@ISA = ("Exporter");
@EXPORT = qw (&changepointCaller);

	sub changepointCaller{

	$SeqFile=$_[0];
	$tabFile=$_[1];
	$Genlen=readGenomeSeq($SeqFile);
	
	open file1, $tabFile
	or die "dead";

	open file2, ">$tabFile.opt"
	or die "dead";

	open file3, ">$tabFile.opt.plot"
	or die "dead";

	$count=`grep -c misc_feature $tabFile`;
	$st=1;	
	$to_prev=0;

	while(<file1>){  								
	chomp($_);
	
		if (m#(\d+)\.\.(\d+)#){
			$left=$1;
			$right=$2;
			$c++;
			print "\n optimizing prediction $c out of $count\n";
		
			#it calls changepoint	
			($newLeft,$newRight)=changepoint($left,$right,$tabFile);
			print " opt boundaries: $newLeft..$newRight\n";
			
			$from=$newLeft;
			$to=$newRight;
			
			#if overlap after the HMM optimization of two consequtive regions
			#then re-update the plot file not to print extra bases
			if($to_prev>=$from){
			$from=$to_prev+1;
			$newLeft=$to_prev+1;
			}
			$to_prev=$to;
		}
				
		if (m#(colour=\d+ \d+ \d+)#){
		$colour=$1;
		}
		if (m#(note=.+)#){
		$note=$1;
		}
		if (m#score=(\d+\.*\d+)#){
		$score=$1;
		
print file2 "FT   misc_feature    $newLeft..$newRight
FT                   /$colour
FT                   /algorithm=\"alien_hunter\"
FT                   /$note
FT                   /score=$score\n";
			
			#begin -> from
			for($i=$st;$i<$from;$i++){
			print file3 "0\n";
			}
			#from -> to
			for($i=$from;$i<=$to;$i++){
			print file3 "$score\n";
			}
			$st=$to+1;
		
		}
		
	}#while
	
		#to -> end		
		for($i=$st;$i<=$Genlen;$i++){
		print file3 "0\n";
		}
	
	close file1;
	close file2;
	close file3;
		
	return (); 

	}
	1;