[Bio] / FigKernelScripts / CSA_first_pass.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/CSA_first_pass.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (view) (download) (as text)

1 : overbeek 1.1 use strict;
2 : overbeek 1.3 #
3 :     # This is a SAS Component
4 :     #
5 : overbeek 1.1 use Data::Dumper;
6 :     use Carp;
7 :     use gjoseqlib;
8 :     use SeedEnv;
9 :    
10 :     my $usage = "usage: CSA_first_pass WorkingDir";
11 :    
12 : overbeek 1.2
13 : overbeek 1.1 my($dir);
14 :     (
15 : overbeek 1.2 ($dir = shift @ARGV)
16 : overbeek 1.1 )
17 :     || die $usage;
18 :    
19 : overbeek 1.2 (-s "$dir/matches") || die "$dir/matches is empty";
20 :    
21 :    
22 : overbeek 1.1 open(MATCHES,"<$dir/matches") || die "could not open $dir/matches";
23 :     my @matches = map { chop;
24 :     my($kmer,$contig1,$strand1,$off1,$contig2,$strand2,$off2) = split(/\t/,$_);
25 :     my $loc1 = [$contig1,
26 :     $off1,
27 :     ($strand1 eq "+") ? ($off1 + (length($kmer)-1)) : ($off1 - (length($kmer)-1))];
28 :     my $loc2 = [$contig2,
29 :     $off2,
30 :     ($strand2 eq "+") ? ($off2 + (length($kmer)-1)) : ($off2 - (length($kmer)-1))];
31 :     [$loc1,$loc2,'pin']
32 :     } <MATCHES>;
33 :     close(MATCHES);
34 :    
35 :     while (my $pin = shift @matches)
36 :     {
37 :     # print STDERR &Dumper(['processing pin',$pin]);
38 :     my $merged;
39 :     while ((@matches > 0) && ($merged = &can_merge($pin,$matches[0])))
40 :     {
41 :     # print STDERR &Dumper(['merged',$matches[0],$pin,$merged]);
42 :     $pin = $merged;
43 :     shift @matches;
44 :     }
45 :     # print STDERR &Dumper($pin);
46 :     print_connection($pin);
47 :     }
48 :    
49 :     sub can_merge {
50 :     my($pin1,$pin2) = @_;
51 :    
52 :     # print STDERR &Dumper(['can_merge',$pin1,$pin2]);
53 :     my($locA1,$locA2,$typeA) = @$pin1;
54 :     my($contigA1,$bA1,$eA1) = @$locA1;
55 :     my($contigA2,$bA2,$eA2) = @$locA2;
56 :    
57 :     my($locB1,$locB2,$typeB) = @$pin2;
58 :     my($contigB1,$bB1,$eB1) = @$locB1;
59 :     my($contigB2,$bB2,$eB2) = @$locB2;
60 :    
61 :     my $merged;
62 :     if ($merged = &can_extend_block($pin1,$pin2))
63 :     {
64 :     return $merged;
65 :     }
66 :     else
67 :     {
68 :     return undef;
69 :     }
70 :     }
71 :    
72 :     sub can_extend_block {
73 :     my($pin1,$pin2) = @_;
74 :    
75 :     my($locA1,$locA2,$typeA) = @$pin1;
76 :     my($contigA1,$bA1,$eA1) = @$locA1;
77 :     my($contigA2,$bA2,$eA2) = @$locA2;
78 :    
79 :     my($locB1,$locB2,$typeB) = @$pin2;
80 :     my($contigB1,$bB1,$eB1) = @$locB1;
81 :     my($contigB2,$bB2,$eB2) = @$locB2;
82 :    
83 :     if (($typeA =~/pin|block/) &&
84 :     ($typeB =~/pin|block/) &&
85 :     ($contigA1 eq $contigB1) &&
86 :     ($contigA2 eq $contigB2) &&
87 :     (&strand($locA1) eq &strand($locB1)) &&
88 :     (&strand($locA2) eq &strand($locB2)))
89 :     {
90 :     my $ext;
91 :     if ((&SeedUtils::between($bA1,$bB1,$eA1) && (! &SeedUtils::between($bA1,$eB1,$eA1))) &&
92 :     (&SeedUtils::between($bA2,$bB2,$eA2) && (! &SeedUtils::between($bA2,$eB2,$eA2))) &&
93 :     (($ext = &SeedUtils::min(abs($eB1-$bA1),abs($eB1-$eA1))) == (&SeedUtils::min(abs($eB2-$bA2),abs($eB2-$eA2)))))
94 :     {
95 :     return [[$contigA1,$bA1,$eB1],[$contigA2,$bA2,$eB2],(($typeA eq "pin") &&
96 :     ($typeB eq "pin") &&
97 :     ($ext == 1)) ? 'pin' : 'block'];
98 :     }
99 :     elsif ((! &SeedUtils::between($bA1,$bB1,$eA1) && (! &SeedUtils::between($bA1,$eB1,$eA1))) &&
100 :     (($ext = ($eB1 - $eA1)) == abs($eB2 - $eA2)) && ($ext < 200))
101 :     {
102 :     return [[$contigA1,$bA1,$eB1],[$contigA2,$bA2,$eB2],'block'];
103 :     }
104 :     }
105 :     return undef;
106 :     }
107 :    
108 :    
109 :     sub strand {
110 :     my($loc) = @_;
111 :    
112 :     return ($loc->[1] <= $loc->[2]) ? '+' : '-';
113 :     }
114 :    
115 :    
116 :     sub print_connection {
117 :     my($x) = @_;
118 :    
119 :     my($loc1,$loc2,$type) = @$x;
120 :     print join("\t",(@$loc1,@$loc2,$type)),"\n";
121 :     }
122 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3