view release on metacpan or  search on metacpan

bin/gbrowse_import_ucsc_db.pl  view on Meta::CPAN

    print STDERR "usage: $0 <UCSC data source name>\n";
    print STDERR "Run $0 --help for details.\n";
    print STDERR "Run $0 --list for list of data sources.\n";
    exit -1;
} else {
    $dbh->do("use $dsn") or die "Could not access $dsn database. Run this script with --list to see valid database names.\n";
}

my $conf_dir = GBrowse::ConfigData->config('conf');
my $data_dir = GBrowse::ConfigData->config('databases');

print STDERR "** During database creation, you may be asked for your password in order to set file permissions correctly.\n\n";

my $dir       = create_database_dir($data_dir,$dsn);
my $scaffolds = create_scaffold_db($dir,$dsn);
my $genes     = create_gene_db($dir,$dsn);
create_conf_file($conf_dir,$dsn,$description,$scaffolds,$genes);
create_source($conf_dir,$dsn,$description);

print STDERR <<END;

** These files have been created for you:

 $scaffolds -- database of chromosome sizes and sequences
 $genes -- database of genes
 $conf_dir/${dsn}.conf -- track configuration file for these databases
 $conf_dir/GBrowse.conf -- updated data source configuration file
END
    ;

if (-x '/usr/sbin/service') {
    print STDERR "\n ** Restarting apache. You may be asked for your password.\n";
    system "sudo service apache2 restart";
} elsif (-x '/etc/init.d/apache2') {
    print STDERR "\n ** Restarting apache. You may be asked for your password.\n";
    system "sudo /etc/init.d/apache2 restart";
} else {
    print STDERR "\n** Please restart apache or other web server.\n";
}

1;

exit 0;

sub print_sources {
    my $dbh = shift;
    my $s = $dbh->selectcol_arrayref('show databases');
    print STDERR join("\n",grep {!/information_schema/} @$s),"\n";
}

sub create_database_dir {
    my ($data_dir,$dsn) = @_;
    my $uid   = $<;
    my ($gid) = $( =~ /^(\d+)/;
    my $dir = "$data_dir/$dsn";
    unless (-d $dir) {
	print STDERR "Creating database directory for $dsn. You may be prompted for your password.\n";
	system "sudo mkdir -p $dir";
    }
    unless (-w $dir) {
	system "sudo chown $uid $dir";
	system "sudo chgrp $gid $dir";
    }
    return $dir;
}

sub create_scaffold_db {
    my ($dir,$dsn) = @_;
    my $path = "$dir/chromosomes";
    mkdir $path unless -e $path;
    open my $db,">$path/chrom_sizes.gff3" or die "$path: $!";
    print $db <<END;
##gff-version 3

END

    print STDERR "Fetching chromosome sizes...\n";

    my $query = $dbh->prepare('select chrom,size from chromInfo order by size')
	or die $dbh->errstr;
    $query->execute;
    my @chroms;
    while (my($chrom,$size) = $query->fetchrow_array) {
	$chrom =~ s/^chr// if $REMOVE_CHR;
	print $db join("\t",
		      $chrom,
		      $dsn,
		      'chromosome',
		      1,
		      $size,
		      '.','.','.',
		      "ID=$chrom;Name=$chrom"),"\n";
	push @chroms,$chrom;
    }

    $query->finish;
    close $db;

    print STDERR "Fetching FASTA files...";
    $ENV{FTP_PASSIVE}=1 unless exists $ENV{FTP_PASSIVE};
    my $prefix = $REMOVE_CHR ? 'chr' : '';
    for my $chr (sort @chroms) {
	my $url  = "ftp://hgdownload.cse.ucsc.edu/goldenPath/$dsn/chromosomes/$prefix$chr.fa.gz";
	my $file = "$path/$prefix$chr.fa.gz";
	print STDERR "$chr...";
	my $code = mirror($url=>$file);
	warn "Fetch of $url returned error code $code\n"
	    if is_error($code);
    }
    print STDERR "done\n";
    print STDERR "Unpacking FASTA files...";
    unlink "$path/chromosomes.fa";
    for my $chr (sort @chroms) {
	my $command = $REMOVE_CHR ? "gunzip -c $path/$prefix$chr.fa.gz | perl -p -e 's/^>chr/>/' >> $path/chromosomes.fa"
	                          : "gunzip -c $path/$prefix$chr.fa.gz >> $path/chromosomes.fa"
	    unless -e "$path/$chr.fa" && -M "$path/chromosomes.fa" > -M "$path/$prefix$chr.fa.gz";
	system $command;
    }
    print STDERR "done\n";

    print STDERR "Creating FASTA index...";
    my $index = Bio::DB::Fasta->new($path) or die "Couldn't create index";
    print "done\n";

    my $wwwuser = GBrowse::ConfigData->config('wwwuser');
    system "sudo chown -R $wwwuser $path";
    return $path;
}

sub create_gene_db {
    my ($dir,$dsn) = @_;
    my $path = "$dir/refGenes";
    mkdir $path unless -e $path;

    my $src_path = "$path/genes.gff3";
    my $db_path = "$path/genes.sqlite";

    open my $db,'>',$src_path or die "$path: $!";
    print $db <<END;
##gff-version 3

END
    print STDERR "Fetching genes...";
    my $query;
    eval {
	$query = $dbh->prepare('select * from refFlat order by geneName') or die $dbh->errstr;
	$query->execute or die $dbh->errorstr;
    } || eval {
	$query = $dbh->prepare('select name2,name,chrom,strand,txStart,txEnd,cdsStart,cdsEnd,exonCount,exonStarts,exonEnds from ensGene order by name2') or die $dbh->errstr;
	$query->execute or die $dbh->errorstr;
    };
    die $@ if $@;
    my $writer = GFFWriter->new($db,$dsn);
    while (my $row = $query->fetchrow_arrayref) {
	$writer->write_transcript($row);
    }
    $query->finish;
    $writer->finish;
    print STDERR "done\n";

    close $db;

    print STDERR "Indexing...";
    system "bp_seqfeature_load.pl -f -c -a DBI::SQLite -d $db_path $src_path";
    print STDERR "done\n";

    return $db_path;
}

sub log10 { log(shift())/log(10)}

sub create_conf_file {
    my ($conf_dir,$dsn,$description,$scaffolds,$genes) = @_;
    my $conf_path = "$conf_dir/${dsn}.conf";
    create_writable_file($conf_path);

    # figure size of chromosomes
    open my $fh1,"$scaffolds/chrom_sizes.gff3" or die "$scaffolds/chrom_sizes.gff3: $!";
    my $max_chrom = 0;
    my ($first_chrom,$first_size);
    while (<$fh1>) {
	my ($chr,undef,undef,undef,$size) = split /\s+/;
	next unless $size;
	$max_chrom = $size if $max_chrom < $size;
	$first_chrom ||= $chr;
	$first_size  ||= $size;

bin/gbrowse_import_ucsc_db.pl  view on Meta::CPAN

glyph        = dna
global feature = 1
database     = scaffolds
height       = 40
do_gc        = 1
gc_window    = auto
strand       = both
fgcolor      = red
axis_color   = blue

[DNA/GC Content:1000000]
hide = 1

[TranslationR]
glyph        = translation
global feature = 1
database     = scaffolds
height       = 20
fgcolor      = blue
strand       = -1
translation  = 3frame
key          = 3-frame translation (reverse)

[TranslationR:1000000]
hide = 1

END
    close $fh3;
}

sub create_source {
    my ($conf_dir,$dsn,$description) = @_;
    my $path = "$conf_dir/GBrowse.conf";
    open my $fh,$path or die "$path: $!";
    my $foundit;
    while (<$fh>) {
	$foundit++ if /\[$dsn\]/;
    }
    close $fh;
    return if $foundit;

    create_writable_file($path);
    open my $fh2,'>>',$path or die "$path: $!";
    print $fh2 "\n";
    print $fh2 <<END;
[$dsn]
description = $description
path        = $conf_dir/${dsn}.conf
END
    ;
    close $fh2;
}

sub create_writable_file {
    my $file = shift;
    return if -e $file && -w $file;

    my $uid   = $<;
    my ($gid) = $( =~ /^(\d+)/;
    system "sudo touch $file";
    system "sudo chown $uid $file";
    system "sudo chgrp $gid $file";
    system "chmod +w $file";
}

package GFFWriter;

sub new {
    my $class = shift;
    my ($fh,$dsn)  = @_;
    return bless {
	fh => $fh,
	dsn => $dsn,
	gene_id        => 'g000000',
	transcript_id  => 't00000',
	last_gene_name => '',
	last_gene_id   => '',
	last_gene      => {},
    },ref $class || $class;
}

# This subroutine is amazingly long and complicated looking, but probably correct
sub write_transcript {
    my $self = shift;
    my $fields = shift;

    my ($gene_name,$accession,$chrom,$strand,$txStart,$txEnd,$cdsStart,$cdsEnd,$exons,$exonStarts,$exonEnds) = @$fields;
    my ($utr5_start,$utr5_end,$utr3_start,$utr3_end,$gid,$tid);
    $gene_name ||= $accession;
    $chrom =~ s/^chr// if $REMOVE_CHR;

    if ($self->{last_gene_name} ne $gene_name
	||
	$self->{last_gene}{chr} ne $chrom                      # avoid some gene name collisions
	||
	abs($self->{last_gene}{start} - $txStart) > 6_000_000  # avoid some gene name collisions
	||
	$self->{last_gene}{strand} ne $strand) 
    {
	$self->write_last_gene;
	$self->{last_gene}           = {};
	$self->{last_gene_name}      = '';
	$gid = $self->{last_gene_id} = $self->{gene_id}++;
    } elsif ($self->{last_gene_id}) {
	$gid = $self->{last_gene_id};
    } else {
	$gid = $self->{last_gene_id} = $self->{gene_id}++;
    }

    $tid = $self->{transcript_id}++;

    my $ORIGIN = 1;
    my $SRC    = $self->{dsn};
    my $fh     = $self->{fh};

    # adjust for Jim's 0-based coordinates
    $txStart++;
    $cdsStart++;

    $txStart  -= $ORIGIN;
    $txEnd    -= $ORIGIN;