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

Annotation of /FigKernelScripts/CSA_get_repeat_dna.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 : overbeek 1.2 #
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_get_repeat_dna WorkingDir";
11 :     $| = 1;
12 :    
13 :     my($dir);
14 :     (
15 :     ($dir = shift @ARGV)
16 :     )
17 :     || die $usage;
18 :    
19 :     my %contigs1 = map { $_->[0] => $_->[2] } &gjoseqlib::read_fasta("$dir/contigs1");
20 :     my %contigs2 = map { $_->[0] => $_->[2] } &gjoseqlib::read_fasta("$dir/contigs2");
21 :    
22 :     system "svr_big_repeats -i 90 -f $dir/contigs1 > $dir/big.repeats1";
23 :     system "svr_big_repeats -i 90 -f $dir/contigs2 > $dir/big.repeats2";
24 :    
25 :     &get_dna_for_repeats("$dir/big.repeats1","$dir/big.repeats1.fasta",\%contigs1);
26 :     &get_dna_for_repeats("$dir/big.repeats2","$dir/big.repeats2.fasta",\%contigs2);
27 :    
28 :     system "formatdb -p F -i $dir/contigs2";
29 :     my @blast_against_ref = grep { ($_ =~ /\s+(\S+)\s+\S+$/) && ($1 < 1.0e-100) }
30 :     `blastall -FF -m 8 -p blastn -d $dir/contigs2 -i $dir/big.repeats1.fasta`;
31 :     open(REP,">$dir/blast.rep1.contigs2") || die "could not open $dir/blast.rep1.contigs2";
32 :     foreach $_ (@blast_against_ref)
33 :     {
34 :     print REP $_;
35 :     }
36 :     close(REP);
37 :    
38 :     sub get_dna_for_repeats {
39 :     my($repeats,$seqF,$contigs) = @_;
40 :    
41 :     open(DNA,">$seqF") || die "could not open $seqF";
42 :     foreach $_ (`cat $repeats`)
43 :     {
44 :     chop;
45 :     my(undef,undef,$contig1,$beg1,$end1,$contig2,$beg2,$end2) = split(/\t/,$_);
46 :     my $seq = &dna_seq($contig1,$beg1,$end1,$contigs);
47 :     $seq || die "$contig1 $beg1 $end1";
48 :     print DNA ">$contig1\_$beg1\_$end1\n$seq\n";
49 :    
50 :     my $seq = &dna_seq($contig2,$beg2,$end2,$contigs);
51 :     $seq || die "$contig2 $beg2 $end2";
52 :     print DNA ">$contig2\_$beg2\_$end2\n$seq\n";
53 :     }
54 :     close(DNA);
55 :     system "formatdb -p F -i $seqF";
56 :     }
57 :    
58 :     sub dna_seq {
59 :     my($contig,$beg,$end,$contigs) = @_;
60 :     if ($beg < $end)
61 :     {
62 :     return substr($contigs->{$contig},$beg-1,$end-($beg-1));
63 :     }
64 :     else
65 :     {
66 :     return &SeedUtils::rev_comp(substr($contigs->{$contig},$end-1,$beg-($end-1)));
67 :     }
68 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3