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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3