get DNA sequence through UCSC DAS server

[[email protected]]

Hi All,

I am trying out the perl code from UCSC help desk posted on biostar to retrieve DNA sequence from UCSC DAS server. After run this script, I didn't get any sequence at all. Could someone have a quick look and point out what might be wrong?

Many thanks!

Stuart
****************************

zzong@bugaboo:~/programming> cat retrieve_dna_from_ucsc.pl
use strict;
use Bio::das; #actually Das

#### USAGE
unless(@ARGV == 2) {
die("USAGE: $0 | Input BED file | Output FASTA file\n\n");
}

#### Access files
open(INPUT, "<$ARGV[0]") or die("Could not open input!\n");
open(OUTPUT, ">$ARGV[1]") or die("Could not open output!\n");

#### DAS db
my $das = Bio::das->new(-source => 'http://genome.cse.ucsc.edu/cgi-bin/das/', -dsn=>'hg19');

#### Line counter
my $count = 0;

#### Work thru each line of BED file
while(defined(my $line = <INPUT>)) {
#### Make sure line begins with chr
unless($line =~ /^chr/) { next }

#### Count lines
$count++;
print "1$count $line\n";

#### Split line
my ($chr, $start, $end, $label) = '';
my @line_bits = ();

@line_bits = split(/\t/, $line);

$chr = $line_bits[0];
$start = $line_bits[1];
$end = $line_bits[2];
$label = $line_bits[3];
chomp($label);

#### Define segment of genome
my $segment = $das->segment("$chr\:$start\,$end");

#### Get DNA
my $dna = $segment->dna;

#### Print sequence to output
print OUTPUT "\>$chr:$start-$end\_$label\n$dna\n";
}

#### Close files
close(INPUT);
close(OUTPUT);

exit;

*******************

[mastal]

get DNA sequence through UCSC DAS server

Did you get any error messages when you ran the script?

[[email protected]]

I figured this out. Minor code changes. The only thing needs to keep in mind is that the input file needs to be tab delimited and every line starts with "chr***.

Regards!

Stuart

**********************

use strict;
use Bio:as;

#### USAGE
unless(@ARGV == 2) {
die("USAGE: $0 Input BED file Output FASTA file\n\n");
}

#### Access files
open(INPUT, "<$ARGV[0]") or die("Could not open input!\n");
open(OUTPUT, ">$ARGV[1]") or die("Could not open output!\n");

#### DAS db

#my $das = Bio:as->new(-source => 'http://genome.cse.ucsc.edu/cgi-bin/das/', -dsn=>'hg19');
my $das = Bio:as->new(-source => 'http://genome.ucsc.edu/cgi-bin/das/', -dsn=>'hg19');
#### Line counter
my $count = 0;

#### Work thru each line of BED file
while(defined(my $line = <INPUT>)) {
#### Make sure line begins with chr
unless($line =~ /^chr/) { next }

#### Count lines
$count++;

#print "1$count $line\n";
print "$count $line\n";
#### Split line
my ($chr, $start, $end, $label) = '';
my @line_bits = ();

@line_bits = split(/\t/, $line);

$chr = $line_bits[0];
$start = $line_bits[1];
$end = $line_bits[2];
$label = $line_bits[3];
chomp($label);

# print "$chr is","$start is","$end is", "$label is";
#### Define segment of genome
my $segment = $das->segment("$chr\:$start\,$end");
#print $segment;
#### Get DNA
my $dna = $segment->dna;
#print $dna;
#### Print sequence to output
print OUTPUT "\>$chr:$start-$end\_$label\n$dna\n";
#print OUTPUT "\>$chr:$start-$end\_$label\n";
}

#### Close files
close(INPUT);
close(OUTPUT);

exit;