ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/chrange.pl
Revision: 24
Committed: Tue Jul 26 21:46:39 2011 UTC (8 years, 1 month ago) by gpertea
File size: 1684 byte(s)
Log Message:
Line File contents
1 #!/usr/bin/perl
2 use strict;
3 use Bio::DB::Fasta;
4 use Getopt::Std;
5 use FindBin;use lib $FindBin::Bin;
6 my $gdefault='/fs/szdata/genomes/human_hg19/all_chrs.fa';
7 my $usage = qq/Usage:
8 chrange.pl [-g <genome_fasta>] [-C][-F][-U] <chromosome>:<coord_start>-<coord_end>
9 Options:
10 -C reverse complement the extracted sequence
11 -U convert sequence to uppercase
12 -F output fasta formatted sequence
13 -g provide the custom genome file instead of the default:
14 $gdefault
15 /;
16 umask 0002;
17 getopts('BCUFg:o:') || die($usage."\n");
18 die($usage."\n") unless @ARGV>0;
19 my $outfile=$Getopt::Std::opt_o;
20 my $gfile=$Getopt::Std::opt_g || $gdefault;
21 die("Error: no such file: $gfile\n") unless -f $gfile;
22 my ($fasta, $upcase, $revcompl) = ($Getopt::Std::opt_F, $Getopt::Std::opt_U, $Getopt::Std::opt_C);
23
24 my %opts = (-maxopen=>1);
25 $opts{-reindex}=1 if $Getopt::Std::opt_B;
26 my $db = Bio::DB::Fasta->new($gfile, %opts)
27 || die("Error creating Bio::DB::Fasta($gfile)\n");
28 #print STDERR " ..fasta file opened.\n";
29 foreach my $gloc (@ARGV) {
30 my ($chr, $start, $rangetype, $stop)=($gloc=~m/^([\w\-\|]+)\:(\d+)([\.\-_\,\:]+)(\d+)/);
31 die($usage."\n") unless $chr && $start>0 && $stop>0;
32 $stop=$start+$stop-1 if ($rangetype eq ':');
33 my $seq=$db->seq($chr.':'.$start.','.$stop);
34 die("No sequence extracted at $chr:$start-$stop\n") unless $seq;
35 $seq=uc($seq) if $upcase;
36 $seq=revCompl($seq) if $revcompl;
37 if ($fasta) {
38 print ">$chr|$start"."_$stop\n".
39 join("\n", unpack('(A70)*', $seq))."\n";
40 }
41 else {
42 print "$seq\n";
43 }
44 } # for each range given
45
46 sub revCompl {
47 my $s=reverse($_[0]);
48 $s =~ tr/AaCcTtGgUuMmRrWwSsYyKkVvHhDdBb/TtGgAaCcAaKkYyWwSsRrMmBbDdHhVv/;
49 return $s;
50 }

Properties

Name Value
svn:executable *