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

Annotation of /FigKernelScripts/CSA_second_pass.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     use Data::Dumper;
3 :     use Carp;
4 :     use gjoseqlib;
5 :     use SeedEnv;
6 :    
7 :     my $usage = "usage: CSA_second_pass WorkingDir";
8 :     $| = 1;
9 :    
10 :     my($dir);
11 :     (
12 :     ($dir = shift @ARGV) && (-s "$dir/matches")
13 :     )
14 :     || die $usage;
15 :    
16 :     my %contigs1 = map { $_->[0] => $_->[2] } &gjoseqlib::read_fasta("$dir/contigs1");
17 :     my %contigs2 = map { $_->[0] => $_->[2] } &gjoseqlib::read_fasta("$dir/contigs2");
18 :    
19 :     my @first_pass = map { chop; [split(/\t/,$_)] } `cat $dir/output.first.pass`;
20 :    
21 :     my $x;
22 :     while ($x = shift @first_pass)
23 :     {
24 :     # print &Dumper(['processing',$x]);
25 :     &print_connection($x);
26 :     my $gap;
27 :     if ((@first_pass > 0) &&
28 :     ($gap = &gen_alignments($x,$first_pass[0],\%contigs1,\%contigs2)))
29 :     {
30 :     # print &Dumper(['gap',$gap]);
31 :     &print_connections($gap);
32 :     }
33 :     }
34 :    
35 :     sub print_connections {
36 :     my($x) = @_;
37 :    
38 :     foreach my $ali (@$x)
39 :     {
40 :     &print_connection($ali);
41 :     }
42 :     }
43 :    
44 :     sub print_connection {
45 :     my($x) = @_;
46 :     print join("\t",@$x),"\n";
47 :     }
48 :    
49 :     sub gen_alignments {
50 :     my($x,$y,$contigs1,$contigs2) = @_;
51 :    
52 :     # print &Dumper(['gen_alignments',$x,$y]);
53 :     my($ln1,$ln2);
54 :     if (($x->[0] eq $y->[0]) && ($x->[3] eq $y->[3]) # if same contigs
55 :     && ($y->[1] > $x->[2]) && (($ln1 = ($y->[1] - ($x->[2]+1))) <= 10000) && ($ln1 > 0))
56 :     {
57 :     if ((((&strand($x->[4],$x->[5]) eq "+") &&
58 :     (&strand($y->[4],$y->[5]) eq "+") &&
59 :     (($ln2 = ($y->[4] - ($x->[5]+1))) > 0)) ||
60 :     ((&strand($x->[4],$x->[5]) eq "-") &&
61 :     (&strand($y->[4],$y->[5]) eq "-") &&
62 :     (($ln2 = ($x->[5] - ($y->[4]+1))) > 0))) &&
63 :     ($ln2 > 0) && (abs($ln2-$ln1) <= (0.1 * $ln1)))
64 :     {
65 :     my $seq1 = &dna_seq($x->[0],$x->[2]+1,$y->[1]-1,$contigs1);
66 :     if (length($seq1) != $ln1) { print &Dumper(['bad seq1',$x,$y,$ln1,$seq1]); die '1' }
67 :     my $incr = ($x->[4] < $x->[5]) ? 1 : -1;
68 :     my $seq2 = &dna_seq($x->[3],$x->[5]+$incr,$y->[4]-$incr,$contigs2);
69 :     if (length($seq2) != $ln2) { print &Dumper(['bad seq2',$ln2,$x,$y,$seq2]); die '2' }
70 :     if (abs(length($seq1) - length($seq2)) > (0.1 * $ln1))
71 :     {
72 :     print &Dumper($x,$y,$seq1,$seq2,$x->[3],$x->[5],$incr,$y->[4],$ln1,$ln2);
73 :     die "bad lengths";
74 :     }
75 :     my @ali = &get_alignments($seq1,$seq2,
76 :     $x->[0],$x->[2]+1,$y->[1]-1,,
77 :     $x->[3],$x->[5]+$incr,$y->[4]-$incr);
78 :     if (@ali > 0) { return \@ali }
79 :     }
80 :     }
81 :     return undef;
82 :     }
83 :    
84 :    
85 :     sub get_alignments {
86 :     my($seq1,$seq2,$contig1,$beg1,$end1,$contig2,$beg2,$end2) = @_;
87 :    
88 :     ($seq1 && $seq2) || confess "bad seqs";
89 :     my @ali;
90 :     my $tmp1 = "tmp1.$$.fasta";
91 :     my $tmp2 = "tmp2.$$.fasta";
92 :     open(TMP,">",$tmp1) || die "could not open $tmp1";
93 :     print TMP ">s1\n$seq1\n>s2\n$seq2\n";
94 :     close(TMP);
95 :     # print "running svr_align_seqs -l -z < $tmp1 > $tmp2\n";
96 :     &run("svr_align_seqs -l -z < $tmp1 > $tmp2");
97 :     # print "look in $tmp1 and $tmp2\n";
98 :     my @ali_seqs = &gjoseqlib::read_fasta($tmp2);
99 :     # unlink($tmp1,$tmp2);
100 :    
101 :     my $ln = length($ali_seqs[0]->[2]);
102 :     my $i=0;
103 :     while ($i < $ln)
104 :     {
105 :     while (($i < $ln) &&
106 :     (substr($ali_seqs[0]->[2],$i,1) eq '-') ||
107 :     (substr($ali_seqs[1]->[2],$i,1) eq '-'))
108 :     {
109 :     $i++;
110 :     }
111 :     my($j,$iden);
112 :     $iden = 0;
113 :     my $c1;
114 :     my $c2;
115 :     for ($j = $i; (($j < $ln) &&
116 :     (($c1 = substr($ali_seqs[0]->[2],$j,1)) ne '-') &&
117 :     (($c2 = substr($ali_seqs[1]->[2],$j,1)) ne '-')); $j++)
118 :     {
119 :     if (uc $c1 eq uc $c2) { $iden++ }
120 :     }
121 :    
122 :     my $ln_piece = $j - $i;
123 :     if (($ln_piece > 5) && (($iden/$ln_piece) >= 0.5))
124 :     {
125 :     push(@ali,[$contig1,$beg1+$i,$beg1+($j-1),
126 :     $contig2,
127 :     ($beg2 < $end2) ? ($beg2+$i,$beg2+($j-1)) : ($beg2-$i,$beg2-($j-1)),
128 :     'ali']);
129 :     }
130 :     $i = $j+1;
131 :     }
132 :     return @ali;
133 :     }
134 :    
135 :     sub run {
136 :     my($cmd) = @_;
137 :    
138 :     my $rc = system($cmd);
139 :     if ($rc)
140 :     {
141 :     die "$rc: $cmd failed";
142 :     }
143 :     }
144 :    
145 :     sub strand {
146 :     my($beg,$end) = @_;
147 :     return ($beg < $end) ? '+' : '-';
148 :     }
149 :    
150 :     sub dna_seq {
151 :     my($contig,$beg,$end,$contigs) = @_;
152 :    
153 :     if ($beg < $end)
154 :     {
155 :     return substr($contigs->{$contig},$beg-1,$end-($beg-1));
156 :     }
157 :     else
158 :     {
159 :     return &SeedUtils::rev_comp(substr($contigs->{$contig},$end-1,$beg-($end-1)));
160 :     }
161 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3