package Bio::GMOD::Bulkfiles;
use strict;
=head1 NAME
Bio::GMOD::Bulkfiles -- produce bulk sequence and feature files
from Chado genome database for public distribution .
=head1 ABOUT Bulkfiles
This generates Fasta, GFF, DNA and other bulk genome annotation files at
ftp://flybase.net/genomes/Drosophila_melanogaster/current/ ..
(and other species)
It works with several FlyBase chado dbs, and with SGDLite chado db
Bulkfiles is mostly self-contained, but uses a few
BioPerl parts plus XML::Simple for configuration files. All of
the organism/database-specific logic should be in these configuration
files (see GMODTools/conf/bulkfiles/)
fbbulk-r4.xml, sgdbulk1.xml .. -- organism/database/release specific options
chadofeatsql.xml -- chado db sql calls to dump features
chadofeatconv.xml -- feature conversion options
=head1 INPUTS
Chado genome databases available (aug 2004) are
ftp://flybase.net/genomes/Drosophila_melanogaster/current/pgsql/chado_r3_2_26_s.gz
http://sgdlite.princeton.edu/download/sgdlite/2004_05_19_sgdlite.sql.gz
Given a postgres db installation, load dump file with
createdb sgdlite_20040519
(zcat 2004_05_19_sgdlite.sql.gz | psql -d sgdlite_20040519 -f - ) >& log.load &
=head1 OUTPUTS
DNA files (full chromosomes) in raw and fasta formats
GFF (v3) feature files
FFF (v1) feature files (used in FlyBase, each complex feature one line
using GenBank/EMBL locations)
Fasta sequence for each selected feature set,
with headers from feature files
BLAST Index files (NCBI)
Chado database overview tables.
A standard genome webpage for access to bulk data.
=head1 REQUIREMENTS and INSTALLATION
Uses a few GMOD, BioPerl, other Perl5 modules, including
Bio::GMOD::Config.pm (and included Config2.pm)
XML::Simple
Program looks for conf/ folder with .xml files.
You likely only need to edit fbbulk-r4.xml equivalent.
=head1 USAGE
# see bin/bulkfiles.pl
use Bio::GMOD::Bulkfiles;
my $bulkfiles= Bio::GMOD::Bulkfiles->new(
configfile => 'fbbulk-r4', # data-release config file, required param
debug => 1, showconfig => 0,
failonerror => 1,
);
$bulkfiles->dumpFeatures(); # extract feature tables from chado sql db
$bulkfiles->sortNSplitByChromosome(); # collate and separate by chromosome
$bulkfiles->dumpChromosomeBases(); # extract chromosome dna files
# produce various output format bulk files from above
my $result= $bulkfiles->makeFiles(
formats => [qw(fff gff fasta blast gnomap)],
chromosomes => [qw(all)] );
print STDERR "Bulkfiles done. result=",$result,"\n";
=head1 WHY Bulkfiles?
(rather than using other middleware layers to chado db - chadoxml,
chadodbi, bioperl, ...)
The general logic is
1. dump all chado db features using simple (and quick) sql,
to common intermediate table files, and chromosome dna to raw files.
The feature info is simple: type, location, name/id, and a few
attributes (db_xrefs,..)
2. postprocess these table files to create the various public use
formats (the time-consuming and configurable part), organized
into per-chromosome files.
Here are some reasons we take this approach:
a. using simple sql to dump all db features to intermediate table
allows easy checks that all features get to bulk files
b. simple sql dump is fast (30 - 60 min for full fly genome),
reliable in getting all mapped features by keeping logic simple
c. process table output in stages - better debugging of steps in
process, and can split processing among computers
c1. the stages are loosely coupled - one can go back, tweek
configurations and get a new output w/o redoing the complete
extraction process.
d. convert one common feature table + dna to several output formats
in one step, or repeatedly as needed.
e. combine features from several chado dbs (flybase now has 3 chado dbs
for d.mel genome features), and add other sources like
flybase cytology features.
f. need fairly complex and data specific configurations - moving
that to config files keeps code reusable.
g. each genome chado database has different policy and choices with
respect to feature, vocabulary and other data. A highly configurable
tool, with data extraction and correction methods that are separate
and tunable is needed to adapt to such variation in genome databases.
=head1 CONFIGURATION
These are essential parts to change for a specific
database release. See e.g., conf/sgdbulk1.xml.
-- include config files allow many of these options
can be packaged separately for reuse among main release configs.
E.g. conf/bulkfiles/fbreleases.xml contains release name, date, databse access
and identified by relid value
SGD Chado DB Lite r1 -- your release title
-- various releases and their chado db name
-- database connection information (release.db overrides this db.name)
-- organism info
scer
Saccharomyces cerevisiae
-- configure chado feature sql dump information
chromosomes -- which sql elements to use
features
matches
analysis
-- configure chado dna/residue sql dump information
-- config. for output files; see conf/bulkfiles/filesets.xml
-- feature sets to make fasta bulk files
-- see conf/bulkfiles/featuresets.xml
gene
mRNA
-- more fasta dump information
MORE CONFIG FILES:
chadofeatsql.xml -- chado db sql calls to dump features
sql for various feature classes, including
chromosomes, gene model, matches, analysis,
synteny, with some helper sql views and postprocessing
scripts.
chadofeatconv.xml -- feature conversion options
this is mostly about how to process specific features
(all SO terms). Tied into Bulkfiles::FeatureWriter code.
=head1 NOTES
genomic sequence file utilities, part2;
parts converted from
flybase/work.local/chado_r3_2_26/soft/chado2flat2.pl
flybase/work.local/chado_r3_2_26/soft/chadosql2flatfeat.pl
Find source here
set cvsd=':pserver:anonymous@cvs.sourceforge.net:/cvsroot/gmod'
cvs -d $cvsd login
cvs -d $cvsd co schema/GMODTools
Backup CVS: set cvsd=':pserver:anonymous@eugenes.org:/bio/cvs'
cvs -d $cvsd co -d GMODTools gmod/schema/GMODTools
=head1 QUICK TEST (Postgres active):
# get soft
cvs -d $cvsd co -d GMODTools schema/GMODTools
# load a genome chado db to Postgres
wget http://http://sgdlite.princeton.edu/download/sgdlite/2004_05_19_sgdlite.sql.gz
createdb sgdlite_20040519
(zcat *sgdlite.sql.gz | psql -d sgdlite_20040519 -f - ) >& log.load
# set GMOD_ROOT to here and run default config
cd GMODTools
env GMOD_ROOT=$PWD perl -I./lib/ bin/bulkfiles.pl sgdbulk1
=head1 UPDATES 2007-Oct, version 1.1
- No chromosome/scaffold/golden_path files. This change is needed to handle
partially assembled genomes with many (1000s to 100,000s) of scaffolds.
Flag no_csomesplit=1 to use this (should become default).
- Gene Ontology association file output, see go_association tags in configurations
formatted as per http://geneontology.org/GO.format.annotation.shtml
- Validate main variables in chado database: ${golden_path}, ${seq_ontology}
This step, on now by default, checks that database contains values set in
configuration for chromosome type, sequence ontology CV name, and any other
critical variables. If failed, db is inspected for real values.
- miscellany bugs cured and configuration updates. tables/overview now again
active.
=head2 no_csomesplit update
this change is needed to handle partially assembled genomes (the common
first draft) with many scaffolds (10,000s of < 1MB to 100,000s of
< 100Kb segments)
The old state assumes a handful of full or nearly full chromosomes implicitly
by producing 1 file/golden_path segment, at least as temporary working
files. This is a problem for any file system, where > 1,000 files/folder
is problematic and > 100,000 can be serious. THis isnt needful design,
just how the program grew (for drosophila genome data with 5 golden_path
segments).
First try: use "sum" as special chromosome name, "all" is already used
specially for final output of public files with all csomes in one file.
Main changes are to these subs here:
See filesets.xml configuration with flags no_csomesplit, perchr and makeall
Features:
sortNSplitByChromosome : produce temporary per-csome feature files (tmp/chadofeat-org-csome)
getFeatTableFiles : collect temp per-csome feat. files for further processing
-- this hands data to downstream BulkWriter packages for final file (GFF,..)
creation. Then the filesets perchr and makeall flags control final file set.
>> need to revise downstream packs (FeatureWriter, ...) to not split into
per-csome files (done before 'makeall' option)
Dna:
-- possibly revise here to use chado sql calls instead of temp dna files?
dumpChromosomeBases: create dna string files from chado feature.residues
see ## for making Unknown 'bag' chromosomes from parts
## $csomeset->{$chr}->{parts}
getBasesFromFiles: return data from dna string files
Note use of Bio::GMOD::Bulkfiles::MyLargePrimarySeq for dna files.
>> first try using only chadodb get residues
failed to create chromosome.fasta (no dna/ files). other fasta looked ok.
See FastaWriter:processFasta : has special raw2Fasta for csomes.,
=head2 Gene Ontology association table
* this is in essense done. See filesets.xml: go_association,
and SQL code in chadofeatsqll.xml,
to use, add to site_defaults, go_association
- this should be a relatively simple subclass for producing
GO association files. Need
(1) Chado SQL to produce basic table, per gene feature, of
linked GO terms via what link table(s)?
-- curated GO terms via feature_cvterm : 52929 in _dmel_r4_3_16?
-- also find IEA GO terms via analysisfeature to alt-species proteins
select fc.cvterm_id, count(distinct fc.feature_id) as n_feature,
(select name from cvterm where cvterm_id = fc.cvterm_id) as term
from feature_cvterm fc group by fc.cvterm_id;
-- add restriction to GO CV tables ...
* something weird here, parameciumDB queries slower than _dmel_r4_3_16,
a larger chado db; maybe we loaded it wrong, with char encoding (UTF-8)
instead of SQL_ASCII
(2) perl code to reformat to GO spec.
=head1 TODO
=head2 Gene Pages XML
- see wfleabase.org/search/ with resultant gene page reports
using simple GeneSummary XML
- output from chado using genepage SQL for all properties/gene feature
=head2 New project configurations
- should add user interface software to prompt for the
few basic required variables (chado db name, golden_path,
choices for output, ...) and
create project configuration.xml from template.xml
=head2 Other useful output files?
- BulkWriter.pm now does basic SQL query to table file,
add some other useful table file outputs? gene IDs x synonyms?
=head2 Remove FFF intermediate for GFF only
- drop the now awkward FFF format for GFF as primary
feature format, use to create Fasta from it.
=head1 SEE ALSO
GMOD::Chado::SeqUtils -- older sequence in/out/check methods for Chado DB
=head1 AUTHOR
D.G. Gilbert, 2004, gilbertd@indiana.edu
=head1 METHODS
=cut
#-----------------
#debug#
#use lib("/bio/argos/common/perl/lib", "/bio/argos/common/system-local/perl/lib");
#use lib("/Users/gilbertd/bio/dev/gmod/schema/GMODTools/lib/");
use POSIX;
use FileHandle;
use File::Basename;
use File::Spec::Functions qw/ catdir catfile /;
use File::Path; ## mkpath
use File::Temp qw/ tempfile tempdir /;
use FindBin qw( $RealBin); #? eval
use DBI;
use Bio::Location::Simple;
use Bio::PrimarySeq; # see MyLargePrimarySeq use
use Bio::GMOD::Bulkfiles::BulkWriter;
# see getWriter() - base class BulkWriter, loads any other subclasses
use Bio::GMOD::Bulkfiles::MyLargePrimarySeq;
use Bio::GMOD::Bulkfiles::MySplitLocation;
# see below require Bio::GMOD::Config2; dont use
our $DEBUG = 0;
my $VERSION = "1.1"; # no_csomesplit update
use vars qw/ @ENV_KEYS @featset @allfeats %mapchr_pattern $fndel /;
# special config keys to put in ENV for Config2 reader
# 0710: add more: ENV_default set : golden_path seq_ontology
# too many rel.variables: relid , release_id == rel ??
BEGIN{@ENV_KEYS= qw( species org date title datadir
rel relfull relid release_url release_id release_date
GFF_source taxon
ftp_url golden_path seq_ontology );}
my $defaultconfigfile="bulkfiles"; # was 'sequtil'
my %dnaseqs=(); #? package global - read only BioseqFile
my @defaultformats= qw(overview fff gff fasta tables blast ); # gnomap
sub new
{
my $that= shift;
my $class= ref($that) || $that;
my %fields = @_;
my $self = \%fields; # config should be one
bless $self, $class;
$self->init();
return $self;
}
sub closeit
{
my ($self)= @_;
my $dbh= delete $self->{dbh};
if ($dbh) { $dbh->disconnect(); }
}
sub DESTROY
{
my $self = shift;
$self->closeit();
}
sub init
{
my $self= shift;
# $self->{tag}= 'Bulkfiles' unless (exists $self->{tag} );
## $self->{outh}= *STDOUT unless ( exists $self->{outh} ); ## DROP THIS ??
$DEBUG= $self->{debug} if defined $self->{debug};
my $cdate= POSIX::strftime("%F %T", localtime( $^T ));
$cdate =~ s/[^\w-]/_/g; ## make sure no space in date !
$self->{date}= $cdate;
$self->{failonerror}= 0 unless defined $self->{failonerror};
$self->{skiponerror}= 1 unless defined $self->{skiponerror};
$self->{ignoredbresidues}= 0 unless defined $self->{ignoredbresidues};
$self->{addids}= 0 unless defined $self->{addids};
$self->{config}={} unless defined $self->{config};
$self->{configfile}= $defaultconfigfile unless defined $self->{configfile};
$self->{verbose}=0 unless defined $self->{verbose};
# no_csomesplit is now a call option, but respect only if >0, replacing config value
$self->readConfig($self->{configfile});
# calls initData()
}
=item readConfig($configfile)
read a configuration file - adds to any loaded configs
-- changed location to conf/bulkfiles/ -- read all of these (no)
but need to locate subdir
=cut
sub readConfig
{
my $self= shift;
$self->{config}= $self->callReadConfig(@_);
$self->initData();
}
# dont dupl this elsewhere e.g. BulkWriter.pm subs
sub callReadConfig
{
my $self= shift;
my ($configfile)= @_;
my $returnConfig= {};
my $replacevariables= $DEBUG; #?? make default, then drop rereadConfig
eval {
my $config2= $self->{config2}; # Config2 object, not hash
unless(ref $config2) {
require Bio::GMOD::Config2;
$self->{config2}= $config2= Bio::GMOD::Config2->new( {
searchpath => [ 'GMODTools/conf/bulkfiles', 'conf/bulkfiles', 'conf', ],
debug => $DEBUG,
read_includes => 1, # process include = 'conf.file'
} );
}
my @showtags= ($self->{verbose}) ? qw(name title date about) : qw(title date);
$config2->setargs( showtags => \@showtags );
# $self->{config}=
$returnConfig= $config2->readConfig( $configfile, {}, {}); ## need empty hashrefs!
## added processing of include="include.conf"
unless ($config2->readConfigOk()) {
warn "** Error with configuration files\n";
die if($self->{failonerror} || $self->{automake});
}
if($replacevariables) { # 0710 update
my @keys= @ENV_KEYS;
# also inspect config for simple keys
foreach my $k (keys %$returnConfig) {
my $v= $$returnConfig{$k};
push(@keys, $k) unless(ref($v));
}
foreach my $k (@keys) {
my $v= $$returnConfig{$k};
$ENV{$k}= $v if(defined $v and not defined $ENV{$k});
}
print STDERR "readConfig:pass2 for variables=@keys\n" if $DEBUG;
# readConfig replace=>1 is bad, replaces 1st/primary values with backup values :(
# new updateVariables will touch only ${variable} set
# recurses thru all of config hash to find '${variable}' matching %ENV keys
$returnConfig= $self->{config2}->updateVariables(
$returnConfig, { Variables => \%ENV, replace =>1 },
);
}
print STDERR $config2->showConfig( $returnConfig, { debug => $DEBUG })
if ($self->{showconfig});
};
if ($@) {
my $cf= $self->{config2}->{filename};
warn "Config2: file=$cf; err: $@";
die if ($self->{failonerror} || $self->{automake});
}
return $returnConfig;
}
sub rereadConfig
{
my $self= shift;
print STDERR "rereadConfig\n" if $DEBUG;
my $config= $self->config;
foreach my $k (@ENV_KEYS) { defined $$config{$k} and $ENV{$k}= $$config{$k}; }
my $newconfig= {};
eval {
## need to replace ${variables} with config values
## so need to change $self->{config} hash, not newconfig
## still need rereadConfig after above updateVariables, as init() changes config
## but could use a 2nd updateVariables here instead ??
$self->{config2}->updateVariables( $config,
{ Variables => \%ENV, replace =>1 }, );
# $newconfig= $self->{config2}->readConfig( $self->{configfile},
# { Variables => \%ENV }, $newconfig, );
# $self->config->{doc} = $newconfig->{doc} if $newconfig->{doc};
};
if ($@) {
my $cf= $self->{config2}->{filename};
warn "Config2: file=$cf; err: $@";
}
# print STDERR "new.doc.content=",$newconfig->{doc}->{content},"\n" if $DEBUG;
}
sub updateConfigVars
{
my $self= shift;
my $config= $self->config;
foreach my $k (@ENV_KEYS) { defined $$config{$k} and $ENV{$k}= $$config{$k}; }
eval { $self->{config2}->updateVariables( $config, { Variables => \%ENV, replace =>1 }, ); };
if ($@) { warn "Config2: file=",$self->{config2}->{filename},"; err: $@"; }
}
sub config { return shift->{config}; }
sub getconfig
{
my $self= shift;
my $cf= $self->{config2};
if ($cf && @_) { return $cf->get(@_); }
}
=item get_filename
$fname= get_filename($org, $chr, $featn, $rel, $format)
make standard output file name "${org}_${chr}_${featn}_${rel}.${format}"
See split_filename; need to parse above out of file name, but
chr, featn, rel can have '_'. Use '-' instead and disallow in parts?
Use '.'?
=cut
sub get_filename
{
my $self= shift;
my( $org, $chr, $featn, $rel, $format)= @_;
unless ( $org ) { $org="noname"; }
if ( $chr ) { $chr="$fndel${chr}"; } else { $chr=''; }
if ( $featn ) { $featn="$fndel${featn}"; } else { $featn=''; }
if ( $rel ) { $rel="$fndel${rel}"; } else { $rel=''; }
unless ( $format ) { $format="txt"; }
my $filename="${org}${chr}${featn}${rel}.${format}";
return $filename;
}
=item split_filename
( $org, $chr, $featn, $rel, $format)= split_filename($filename)
parse standard output file name "${org}-${chr}-${featn}-${rel}.${format}"
THIS IS problematic - both featn and rel can have _, need other split method
now chr can have _: dpse_4_group2_transposable_element_rel_2_1..
BAD NEWS: dpse 4_group5; XR_group3; ... underscore wont do it
$fndel=".";
dpse.4_group2.transposable_element.rel_2.1.fasta
$fndel="-";
dpse-4_group2-transposable_element-rel_2.1.fasta
$fndel="=";
dpse=4_group2=transposable_element=rel_2.1.fasta
$fndel="--";
dpse--4_group2--transposable_element--rel_2.1.fasta
$fndel="#";
dpse#4_group2#transposable_element#rel_2.1.fasta
=cut
sub split_filename
{
my $self= shift;
my ($fname,$no_orgchr)= @_;
my( $org, $chr, $featn, $rel, $format, $path, $featORrel, $gz, $xtra);
if ($fname =~ s/(\.gz)$//) { $gz=$1; }
($fname, $path, $format) = File::Basename::fileparse($fname, '\.[^\.]+');
$format .= $gz if $gz; #??
if($no_orgchr) {
# only fname, format valid ...
($org,$chr,$featn,$rel)= ('')x4;
$featn= $fname; #???
$chr="all";
} else {
#my $nu= ($fname =~ tr/$fndel/$fndel/);## bad
my @v= split(/$fndel/, $fname, 4); #! rel can have _; !! chr can have _; feat might have _
my $nu= scalar(@v);
if (@v == 4) { ($org,$chr,$featn,$rel)= @v; $nu++ if (index($rel,$fndel)>=0); }
elsif (@v == 3) {
($org,$chr,$featORrel)= @v; # might also be $feat,$chr,$xxx
if ($featORrel =~ /\d/) { $rel= $featORrel; $featn=''; }
else { $featn= $featORrel; $rel=''; }
}
elsif (@v == 2) { ($org,$chr)= @v; $featn= $rel='';} ## might also be ($featn,$chr)
if ($nu > 4 || $nu<3) {
warn "split_filename : ambiguous parts $org,$chr,$featn,$rel from '$fname'\n";
}
}
return ( $org, $chr, $featn, $rel, $format);
}
=item dnafile($chr)
locate file of raw chromosome dna bases
?? change naming of dnafile() to same as others using get_filename
now: dna-chrIX.fasta,raw
to : scer_chrIII_dna_sgdr1.fasta,raw
=cut
sub dnafile
{
my ($self, $chrOrFile)= @_;
#? search for it .. if ($self->config->{dnafiles}->{path});
if ($chrOrFile) {
my $dnafile="";
if (-e $chrOrFile) { $dnafile= $chrOrFile; }
else {
my $org= $self->config->{org};
my $rel= $self->config->{rel};
my $fn= $self->get_filename($org,$chrOrFile,'dna',$rel,'raw');
# my $fn= "dna-$chrOrFile.raw";
$dnafile = catfile( $self->{dnadir}, $fn);
}
$self->{dnafile}= $dnafile;
}
return $self->{dnafile};
}
=item remapArm
-- unordered contigs -- singles (? no feats) and doubles - put into common out files?
-- if so, need to offset start/end to fit into unorderd 'chromosome'
See also FeatureWriter remapName, remapArm ...
=cut
sub remapArm
{
my ($self,$arm,$fmin,$fmax,$strand,$csomeset)= @_;
my $savearm= $arm;
$csomeset= {} unless($csomeset);
if ($csomeset->{'_part2chr'}->{$arm}) {
$arm= $csomeset->{'_part2chr'}->{$arm};
}
else {
foreach my $mp (sort keys %mapchr_pattern) {
next if ($mp eq 'null'); # dummy?
my $from= $mapchr_pattern{$mp}->{from}; next unless($from);
my $to = $mapchr_pattern{$mp}->{to};
#$arm =~ s/$from/$to/g;
if ($to =~ /\$/) { $arm =~ s/$from/eval($to)/e; }
else { $arm =~ s/$from/$to/g; }
}
}
## drop arm if $arm is mapped to ""
## need to trap $arm not in $csomeset -- errors like Contig_Contig below
## which either need to be remapped or otherwise handled.
##
if (!$csomeset->{$arm}) {
}
elsif ($csomeset->{$arm} && ref($csomeset->{$arm}->{parts})) {
my $parts= $csomeset->{$arm}->{parts};
MATCHPART:
foreach my $p (@$parts) {
my $nm = $p->{name}; #? or id ?
if ($nm eq $savearm) { #?? is this ok
my $b = $p->{start};
my $e = $p->{length} + $b;
my $st= $p->{strand};
if ($st < 0) { #($st eq '-')?? do we need to flip all - min,max relative to arm.e ?
$strand= -$strand;
($fmax,$fmin) = ($e - $fmin-1, $e - $fmax-1);
}
else {
$fmin += $b - 1;
$fmax += $b - 1;
}
last MATCHPART;
}
}
}
return($arm,$fmin,$fmax,$strand,$savearm);
}
=item dumpChromosomeBases
$sequtil->dumpChromosomeBases( \@chromosomes or $config->chromosomes)
foreach chr @chromosomes
write dnafile() getBasesFromDb($chrID)
>> moved out >>$sequtil->raw2Fasta() if $config->{dofasta}; -- write from db to files
oct04 -- add option to create 'chr-U' from bag of Unknown_* golden_path entries
in getChromosomeTable.
** See FeatureWriter remapArm and precursor work w/ Dpse; need to create
ordered Uknown chr, joining contigs end-to-end.
oct07 -- need no_csomesplit / no dna files option (for 100,000s of scaffolds cases)
=cut
sub dumpChromosomeBases
{
my ($self, $chromosomes)= @_;
#? or chromosome->{}
my @files=();
my $csomeset= $self->getChromosomeTable();
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
my $no_csomesplit= $self->config->{no_csomesplit} || 0; # FIXME: 0710
#? $chromosomes=["sum"] if($no_csomesplit); # 2007oct; not "all", other special csome name !
# ^^ can we use one dna.raw file with simple csome byte index instead of many files?
# or chado db feature.residues fetch as csome 'file'
if($no_csomesplit) {
my $chr="sum"; #? "sum" or "all"
my $dnafile= $self->dnafile($chr); #? add $org option
return [ { path => $dnafile, type => 'dna/raw', name => $chr, chr => $chr, } ];
}
##? change naming of dnafile() to same as other files?
## using get_filename
my $saveorg= $self->config->{org};
foreach my $chr (@$chromosomes) {
next if $chr eq 'all';
my $spp= $csomeset->{$chr}->{species}; ## FIX $dnafile name for this !
my $org= $self->speciesAbbrev($spp);
## next if($org && $org ne $self->config->{org});
## skip if $spp ne $self->config->{species} ? or Rename ; but need spp abbrev
$self->config->{org}= $org if $org;
my $dnafile= $self->dnafile($chr); #? add $org option
$self->config->{org}= $saveorg;
print STDERR "dumpChromosomeBases $dnafile\n" if $DEBUG;
if (-e $dnafile) {
warn "dumpChromosomeBases: wont overwrite $dnafile";
# die if failonerror ?? optionally clean/rewrite ?
}
## for making Unknown 'bag' chromosomes from parts
elsif (ref($csomeset->{$chr}->{parts})) {
my $parts= $csomeset->{$chr}->{parts};
my $len= 0; my $np= 0;
open(DNA,">$dnafile");
foreach my $p (@$parts) {
my $id= $p->{id};
my $bases= ($id) ? $self->getBasesFromDb($id,1) : '';
print DNA $bases if ($bases);
$len += length($bases); $np++;
}
close(DNA);
print STDERR " dumped $chr parts=$np, total_length=",$len,"\n" if $DEBUG;
push(@files, { path => $dnafile, type => 'dna/raw', name => $chr, chr => $chr, });
}
else {
my $id= $csomeset->{$chr}->{id} || $chr;
my $bases= $self->getBasesFromDb($id,1);
if ($bases) {
open(DNA,">$dnafile"); print DNA $bases; close(DNA);
print STDERR " dumped length=",length($bases),"\n" if $DEBUG;
push(@files, { path => $dnafile, type => 'dna/raw', name => $chr, chr => $chr, });
}
else {
warn "dumpChromosomeBases: no bases for $dnafile\n";
die if $self->{failonerror};
}
}
}
return \@files;
}
=item getChromosomeFiles()
return fileset of dna/raw chromosomes
=cut
sub getChromosomeFiles
{
my ($self, $chromosomes)= @_;
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
my @files=();
## FIXME: no_csomesplit
my $no_csomesplit= $self->config->{no_csomesplit} || 0; # FIXME: 0710
foreach my $chr (@$chromosomes) {
next if $chr eq 'all';
my $dnafile= $self->dnafile($chr);
if (-e $dnafile) {
push(@files, { path => $dnafile, type => 'dna/raw', name => $chr, chr => $chr });
}
}
return \@files;
}
=item getDumpFiles($targets)
return list of feature dump files
=cut
sub getDumpFiles
{
my ($self, $targets, $fdump)= @_;
$fdump= $self->config->{featdump} unless($fdump);
my @files=();
my @missing=();
my $seqsql = $self->getSeqSql($fdump->{config},$fdump->{ENV});
my $outpath= $self->getReleaseSubdir( $fdump->{'path'} || catdir( $self->config->{TMP}, "featdump") );
$fdump->{'path'}= $outpath; # save for reuse
my $sqltag = $fdump->{tag} || "feature_sql";
my $sqltype = $fdump->{type};
unless($targets) { $targets = $fdump->{target}; } # should be array ?
unless($targets) { my @tg= sort keys %{$seqsql->{$sqltag}}; $targets= \@tg; }
unless(ref $targets) { $targets= [ $targets ]; }
foreach my $tgname (@$targets)
{
my $fs= $seqsql->{$sqltag}->{$tgname};
unless($fs) { next; }
my $type= $fs->{type};
my $sql = $fs->{sql};
my $outn= $fs->{output} || $tgname.".tsv";
unless($sql && (!$sqltype || $type =~ m/\b$sqltype\b/)) {
warn "getDumpFiles skip: $tgname/$type<>$sqltype\n";
next; } #??
my $outf= catfile($outpath,$outn);
if(! -e $outf && -e "$outf.gz") { $outf .=".gz";}
# changed keys: name => to target => ; file => to name =>
if (-e $outf) {
push(@files, { path => $outf, type => $type, target => $tgname, name => $outn, });
}
else {
push(@missing, { path => $outf, type => $type, target => $tgname, name => $outn, });
}
}
return (wantarray) ? (\@files,\@missing) : \@files;
}
=item getFastaFiles()
return fileset of available features/fasta
=cut
sub getFastaFiles
{
my ($self, $chromosomes)= @_;
my @files=();
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
#? drop this; moved fasta out of dnadir ..
foreach my $chr (@$chromosomes) {
my $dnafile= $self->dnafile($chr);
(my $fastafile = $dnafile.".fasta") =~ s/\.raw//;
if(! -e $fastafile && -e "$fastafile.gz") { $fastafile .=".gz";}
if (-e $fastafile) {
push(@files, { path => $fastafile, type => 'chromosome/fasta', name => $chr, chr => $chr});
}
}
my $fset= $self->getFilesetInfo('fasta');
my $fadir= $self->getReleaseSubdir( $fset->{path} || 'fasta/');
if (opendir(D, $fadir)) {
foreach my $fa (grep(/^\w/,readdir(D))) {
my ( $org, $chr, $featn, $rel, $format)= $self->split_filename($fa);
$featn = 'feature' unless($featn);
next unless( grep {$chr eq $_} @$chromosomes ); #?
push(@files,
{ path => "$fadir/$fa", type => "$featn/fasta", name => $fa, chr => $chr });
}
closedir(D);
}
return \@files;
}
sub getFeatFiles
{
my ($self, $chromosomes)= @_;
my @files=();
my $fset= $self->getFilesetInfo('fff');
my $featdir= $self->getReleaseSubdir( $fset->{path} || 'fff/');
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
if (opendir(D, $featdir)) {
foreach my $fa (grep(/^\w/,readdir(D))) {
my ( $org, $chr, $featn, $rel, $format)= $self->split_filename($fa);
next unless( grep {$chr eq $_} @$chromosomes );
$featn ||= 'feature';
push(@files,
{ path => "$featdir/$fa", type => "$featn/fff", name => $fa, chr => $chr });
}
closedir(D);
}
return \@files;
}
sub getFeatTableFiles
{
my ($self, $chromosomes)= @_;
my @files=();
my $fdump= $self->config->{featdump};
my $outpath= $self->getReleaseSubdir( $fdump->{'path'} || "tmp/") ;
my $outname= catfile( $outpath, $fdump->{splitname} || "chadofeat");
my $no_csomesplit= $self->config->{no_csomesplit} || 0; # FIXME: 0710
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
##my $spp= $csomeset->{$chr}->{species}; ## FIX $dnafile name for this !
my $spp= $self->config->{species} || $self->config->{org};
my $orgabbr= $self->speciesAbbrev($spp, "org4");
print STDERR "# feature_table org=$orgabbr species=$spp path=$outpath \n" if $DEBUG;
$chromosomes=["sum"] if($no_csomesplit); # 2007oct; not "all", other special csome name !
foreach my $chr (@$chromosomes) {
next if $chr eq 'all';
my $fn;
#^^ FIXME for species file name: chadofeat-dmel2L.tsv .. chadofeat-dpseXR_group9.tsv
# ? add spp to getChromosomes() ?
$fn= "$outname-$orgabbr$chr.tsv";
if(! -e $fn && -e "$fn.gz") { $fn .=".gz"; }
#$fn= "$outname-$chr.tsv" unless( -e $fn);
push(@files, { path => $fn, type => 'feature_table', name => $chr, org => $orgabbr });
}
print STDERR "# feature_table files=",join(",",map{ basename($_->{path})}@files),"\n" if $DEBUG;
return \@files;
}
sub getFilesetInfo
{
my ($self, $type)= @_;
## regularize configs so new format can be added w/o special tag names
my $fset;
$fset= $self->config->{fileset_override}->{$type};
return $fset if (ref $fset);
$fset= $self->config->{fileset}->{$type};
return $fset if (ref $fset);
my @oldset= qw( featdump dnafiles featfiles fastafiles blastfiles gnomapfiles gbrowsefiles );
foreach my $ms (@oldset) {
my $fset= $self->config->{$ms};
if (ref $fset && $fset->{type} eq $type) { return $fset; }
}
return {};
}
sub getFilesByType
{
my ($self, $type, $chromosomes)= @_;
my @files=();
$chromosomes= $self->getChromosomes() unless (ref $chromosomes);
## $type maybe [] arrayref
my @types= (ref $type) ? @$type : ( $type ) ;
foreach $type (@types) {
my $fset= $self->getFilesetInfo($type);
next unless $fset;
my $path= $fset->{path} || $type;
my $dir = $self->getReleaseSubdir( $path, 'nocreate' );
my $no_orgchr= (defined $fset->{no_orgchr}) ? $fset->{no_orgchr} : 0;
if (opendir(D, $dir)) {
my $filepattern= '\w';
# ($filepattern, undef) = File::Basename::fileparse($path) if ($path !~ m,/$,);
# >> err fileparse(): need a valid pathname
if ($path !~ m,/$,){
my $e= rindex($path,'/');
$filepattern= ($e<0) ? $path : substr($path,$e+1);
}
foreach my $fa (grep(/^$filepattern/,readdir(D))) {
my ( $org, $chr, $featn, $rel, $format);
( $org, $chr, $featn, $rel, $format)=$self->split_filename($fa,$no_orgchr);
next unless( $no_orgchr || grep {$chr eq $_} @$chromosomes );
$featn ||= 'feature';
push( @files,
{ path => "$dir/$fa", type => "$featn/$type",
name => $fa, format => $format,
chr => $chr, rel => $rel, org => $org,
no_orgchr => $no_orgchr,
});
}
closedir(D);
}
}
return \@files;
}
sub getFiles
{
my ($self, $type, $chromosomes)= @_;
my @filesets=();
## $type maybe [] arrayref
my @types= (ref $type) ? @$type : ( $type ) ;
foreach $type (@types) {
my $fileset;
## old way
CASE: {
$type eq 'feature_table' && do { $fileset= $self->getFeatTableFiles($chromosomes); last CASE };
$type eq 'dna' && do { $fileset= $self->getChromosomeFiles($chromosomes); last CASE };
#$type eq 'fff' && do { $fileset= $self->getFeatFiles($chromosomes); last CASE };
#$type eq 'fasta' && do { $fileset= $self->getFastaFiles($chromosomes); last CASE };
# new way:
$fileset= $self->getFilesByType($type, $chromosomes);
}
push(@filesets, @$fileset) if $fileset;
}
return \@filesets;
}
sub gzipFiles
{
my ($self, $formats, $chromosomes)= @_;
my @formats= ref $formats ? @$formats : qw( fff gff fasta dna gnomap blast) ;
foreach my $type (@formats) {
my $finfo= $self->getFilesetInfo($type);
if (ref $finfo && $finfo->{dogzip}) {
my $fileset= $self->getFiles($type, $chromosomes);
print STDERR "gzipping $finfo->{path}\n" if $DEBUG;
foreach my $fs (@$fileset) {
system("gzip -f ".$fs->{path}) if (-e $fs->{path} && $fs->{path} !~ /\.gz$/);
}
}
}
}
#=============================================
=item sortNSplitByChromosome($fileset)
sort chado feature dump fileset by arm, location
and split into chromosome file set
=cut
sub sortNSplitByChromosome
{
my ($self, $fileset)= @_;
my $sorter=`which sort`; chomp($sorter); ## '/bin/sort'; '/usr/bin/sort';
## WATCH OUT - TAB here in ' ' ; does shell understand '^I' instead ?
#DEBUG off#my $sortfeaturescmd= "$sorter -t' ' -k 1,1 -k 2,2n"; #? add -k 3,3rn ; nope end not there
my $sortfeaturescmd= "$sorter -t' ' -k 1,1 -k 2,2n -k 3,3rn";
$fileset= $self->getDumpFiles() unless(ref $fileset);
# return undef unless(ref $fileset);
my $fdump = $self->config->{featdump};
my $outpath= $self->getReleaseSubdir( $fdump->{'path'} || "tmp/") ;
my $outname= catfile( $outpath, $fdump->{splitname} || "chadofeat");
my $intype = $fdump->{type};
my $no_csomesplit= $self->config->{no_csomesplit} || 0; # FIXME: 0710
my $tfset = $self->getFilesetInfo('tables');
my $tabdir = $self->getReleaseSubdir( $tfset->{path} || 'tables/');
my $sumfile= catfile( $tabdir, "feature_table-summary.txt");
## check first existance of outname files, and age
## if newer than input fileset, leave as is, return file names ?
my $chromosomes= $self->getChromosomes(); ## $self->config->{chromosomes};
my $chr= $$chromosomes[0];
##FIXME: my $testout= $outname-$orgabbr.$chr; ## dmel2R,
my $orgabbr= $self->config->{org} || $self->speciesAbbrev("","org4");
my $testout= "$outname-$orgabbr$chr.tsv";
if(! -e $testout && -e "$testout.gz") { $testout .=".gz"; }
my $uptodate= (-e $testout) ? 1 : 0;
my $scmd="";
foreach my $fs (@$fileset) {
my $fp= $fs->{path};
next unless($fs->{type} eq $intype); #??
## if(! -e $fp && -e "$fp.gz") { $fp .=".gz"; } ## fileset will already have .gz if so
unless(-e $fp) { warn "missing dumpfile $fp"; next; } # die if $self->{failonerror};
if ($uptodate && _isold($fp, $testout)) { $uptodate= 0 ; }
$scmd .= "$fp ";
}
if($uptodate) {
my @files=();
foreach my $chr (@$chromosomes) {
my $fn= "$outname-$orgabbr$chr.tsv"; ## dmel2R,
push(@files, { path => $fn, type => 'feature_table', name => $chr, chr => $chr, });
}
return \@files;
}
unless($scmd) { warn "sortNSplitByChromosome: no dumpfiles at $outpath"; return undef; }
# die if $self->{failonerror};
if($scmd =~ /\.gz /){ $scmd= "gunzip -c $scmd"; } else { $scmd= "cat $scmd"; }
$scmd = "$scmd | $sortfeaturescmd |";
print STDERR "sortNSplitByChromosome:\n $scmd\n" if $DEBUG;
print STDERR " to csomeSplit($outname)\n" if $DEBUG;
open(FS,$scmd) || die $scmd;
my $files= $self->csomeSplit(*FS, $outname, $sumfile, $no_csomesplit);
close(FS);
## with remapArm, need to resort each output file
foreach my $finfo (@$files) {
next unless ($finfo->{dosort});
my $path= $finfo->{path};
my $cmd="cat $path | $sortfeaturescmd > $path.new";
if (system($cmd) == 0) { rename("$path.new",$path); $finfo->{dosort}= 0; }
else { my $err= $? >> 8; warn "resort $path failed: $err"; }
}
return $files;
}
=item csomeSplit($inh, $outname, $sumfile)
split sorted input feature_table files into per-chromosome (and now per-organism)
feature fileset
2007Oct: add this flag ENV_default: output_by_golden_path=1 instead use: csome_split=1
to separate output by chromosome/scaffold/..., otherwise lump output by feature types,
for all csomes.
-- how much work to implement this? as all downstream outputs are csome-split now.
=cut
sub csomeSplit
{
my($self, $inh, $outname, $sumfile, $no_csomesplit)= @_;
$outname ||= "chadofeats";
my @files=();
my $fh= undef;
my %fhs=();
my %csomefeats= ();
my $lastoid='';
my $csomeset = $self->getChromosomeTable(); ## pass value
my $organisms= $self->config->{organism};
my %fhsorts=();
while(<$inh>) {
next unless(/^\w/);
chomp();
my ($arm,$fmin,$fmax,$strand,$orgid,$type,$name,$id,$oid,$attr_type,$attribute)
= split("\t");
my $oldarm;
my $oldfmax= $fmax;
# bug 0710: getting 's.cerevisiae' instead of 'scer'; from chado db ??
my $orgabbr=
$organisms->{$orgid}->{org4} ||
$organisms->{$orgid}->{abbreviation};
$orgabbr = "org$orgid" unless($orgabbr); #?? skip null orgabbr : not wanted ?
($arm,$fmin,$fmax,$strand,$oldarm)= $self->remapArm($arm,$fmin,$fmax,$strand,$csomeset); # for Unknown.. and other fragments ? need to do before sorter call
## -- BUG: remapArm U isnt sorted; need to check after 1st create files
## -- BUG2: renameArm needs to be called other places, change ChromosomeTable
next unless($arm && $arm ne 'skip'); #?
if ($attr_type eq 'to_species') {
my $toorg= $self->speciesAbbrev($attribute); ## , "org4"
## my $toorg= $organisms->{$attribute}->{abbreviation};
$attribute= $toorg if $toorg;
}
my $dosort= ($arm ne $oldarm || $fmax ne $oldfmax);
my $armfile= ($no_csomesplit) ? "sum" : $arm;
my $fhname= $orgabbr.$armfile; ## dmel2R,
unless($fhs{$fhname}) {
my $fn= "$outname-$fhname.tsv";
$fh= $fhs{$fhname}= new FileHandle(">$fn");
die "Failed to create csome feature file '$fn'" unless($fh);
push(@files, { path => $fn, name => $armfile, chr => $arm, org => $orgabbr,
type => 'feature_table', # should be $fs->{type} == feature_table
dosort => $dosort,
});
}
$fhsorts{$fhname} += $dosort;
$fh= $fhs{$fhname};
## drop $orgid from this output set
print $fh join("\t",$arm,$fmin,$fmax,$strand,$type,$name,$id,$oid,$attr_type,$attribute)
,"\n";
unless($oid eq $lastoid) {
$csomefeats{$fhname}{$type}++; $csomefeats{all}{$type}++;
}
$lastoid= $oid;
}
foreach my $fhname (keys %fhs) { $fh= $fhs{$fhname}; close($fh) if $fh; }
foreach my $f (@files) {
my $fhname= $f->{org} . $f->{chr};
$f->{dosort}= 1 if ($fhsorts{$fhname});
}
## these counts are bad; include dup rows/feature (e.g 4x gene count)
## need to use distinct oid .
if ( $sumfile ) {
$fh= new FileHandle(">$sumfile");
my $title = $self->config->{title};
my $date = $self->config->{date};
my $org = $self->config->{species} || $self->config->{org};
#^^ use $orgid !??
print $fh "# Database feature summary for $org from $title [$date]\n";
my @fl= grep { 'all' ne $_ } sort keys %csomefeats;
foreach my $arm ('all', @fl) {
print $fh (($arm eq 'all') ? "\n# ALL chromosomes\n" : "\n# Chromosome $arm\n");
foreach my $t (sort keys %{$csomefeats{$arm}}) {
print $fh "$t\t$csomefeats{$arm}{$t}\n";
}
print $fh "#","="x50,"\n";
}
close($fh);
push(@files, { path => $sumfile, type => 'feature/summary', name => 'summary', });
}
return \@files;
}
sub preMake
{
my $self= shift;
my ($what)= @_;
my $ok= 0;
my ($dumpfiles,$chrfeats,$dnafiles);
# parts from bulkfiles.pl not run by default; save user some hassle here
if ($what =~ /feature_table/ && !$self->{didmake}{feature_table}) {
warn "Automaking feature_table files\n";
$dumpfiles= $self->dumpFeatures();
$chrfeats= $self->sortNSplitByChromosome();
$self->{didmake}{feature_table}=1;
$ok= 1 if(@$chrfeats);
}
if ($what =~ /dna/ && !$self->{didmake}{dna}) {
warn "Automaking dna files\n";
$dnafiles= $self->dumpChromosomeBases();
$self->{didmake}{dna}=1;
$ok= 1 if(@$dnafiles);
}
return $ok;
}
sub missingData
{
my $self= shift;
my ($fileset, $maketype, $what)= @_;
if( @$fileset && $self->{automake} ) {
my @didmake= keys %{$self->{didmake}};
my $current= grep /$maketype/, @didmake;
unless($current) {
(my $msg = $what) =~ s/;.*$//;
warn "Using pre-existing $msg\n";
warn "Use '/bin/rm -rf ",$self->getReleaseDir(),"' for clean make\n"
if($what =~ /chromosomes/);
}
}
unless(@$fileset){ warn "Missing ",$what,"\n"; die if $self->{failonerror}; }
}
=item makeFiles( %args )
primary method;
mostly it calls FeatureWriter() package to handle
also prints any config->doc entries
arguments:
infiles => \@fileset, # required
formats => [ 'gff', 'fasta', 'fff', 'blast' ] # optional
=cut
sub makeFiles
{
my $self= shift;
my %args= @_;
my @results=();
print STDERR "makeFiles\n" if $DEBUG; # debug
$self->rereadConfig(); # == updateConfigVars; replace doc ${values}
$self->writeDocs( $self->config->{doc} ); #? unless already wrote ? move this to Writer module
my $automake= $args{automake} || $self->{automake} ;
## dang change in caller bulkfiles.pl; dropped @outformats, now find
## downstream writers look at this input %args for that: should it (a) ignore $args{outformats}
## or (b) should this method add \@outformats to args?
my @outformats=(); # check config
if ($args{formats}) {
my $formats= $args{formats};
@outformats= (ref $formats) ? @$formats : ($formats);
}
unless(@outformats>0) {
@outformats= @{ $self->config->{outformats} || \@defaultformats } ;
}
# automake, etc, add outformat fff if need be
# push(@outformats,"fff") # still a dependency, but intermediate format
# if ($automake && !grep(/fff/, @outformats) && grep(/fasta/, @outformats));
# bug^ recreates fff already made;
$args{formats}= \@outformats; #??? want this
my %outformats= map{ $_,1; } @outformats;
print STDERR "makeFiles: outformats= @outformats\n" if $DEBUG;
if (delete $outformats{'overview'}) {
my $overviewset = $self->getFilesetInfo('overview');
if($overviewset) {
## check for already done overview files ..
my ($ovfiles,$ovmissing)= $self->getDumpFiles(['overview'], $overviewset);
## FIXME: bug in above x config? 'overview' files are tagged 'summary' in chadofeatsql ??
if(1 or scalar(@$ovmissing)) {
$ovfiles = $self->dumpFeatures($overviewset, undef, "colnames");
}
my $ovlist= join(" ",map {$_->{name}} @$ovfiles);
push @results, "overviews:$ovlist";
warn "** Please review overview tables for validity **\noverviews:$ovlist\n" if($automake||$DEBUG); # && return ??
$self->getOrganismTable() if($ovlist =~ /organism/);
}
}
## 0710: insert here optional validateVariables: seq_ontology, golden_path at least
## valid=0 default flag handles this for newuser: $self->config->{newuser}
$self->validateVariables() unless( $self->config->{valid} );
## getChromosomes needs chromosomes.tsv .. essential it exist here
my $fileset= $self->getDumpFiles(['chromosomes']);
if($automake && ! @$fileset) {
$self->preMake('feature_table'); # returns $ok
$fileset= $self->getDumpFiles(['chromosomes']);
}
$self->missingData( $fileset, 'feature_table', "chromosomes.tsv table file; make with -featdump");
my $chromosomes= undef;
if (ref $args{chromosomes}) { $chromosomes= $args{chromosomes}; $args{noall}=1; }
elsif (ref $args{chr}) { $chromosomes= $args{chr}; $args{noall}=1; }
## dont do 'all' if subset !
unless (ref $chromosomes && @$chromosomes > 0 && $chromosomes->[0] ne 'all') {
$chromosomes= [ 'all', @{$self->getChromosomes()} ];
$args{noall}=0;
}
# FIXME trick - getFeatureWriter loads common config for featmap/featset needed by others
my $featwriter= $self->getWriter('fff');
my $featfiles = $self->getFiles('fff', $chromosomes);
if ($automake && !@$featfiles && ($outformats{'fasta'} || $outformats{'gnomap'})) { #(grep /fasta|gnomap/, @outformats))
$outformats{'fff'}=1;
}
## this one takes a while; split chromosomes among processors
if ($outformats{'fff'} || $outformats{'gff'}) { ## grep /fff|gff/, @outformats)
delete $outformats{'fff'}; delete $outformats{'gff'};
my $chrfeats = $self->getFiles('feature_table', $chromosomes);
if ($automake && !@$chrfeats) {
$self->preMake('feature_table'); # returns $ok
$chrfeats = $self->getFiles('feature_table', $chromosomes);
}
if ($DEBUG) { print STDERR "read feature tables= ",join(" ",map {$_->{name}} @$chrfeats),"\n"; }
$self->missingData( $chrfeats, 'feature_table', "feature_table files; make with -featdump");
push @results, $featwriter->makeFiles( %args,
infiles => $chrfeats, chromosomes => $chromosomes );
}
my $featfiles = $self->getFiles('fff', $chromosomes);
my $dnafiles = $self->getFiles('dna', $chromosomes);
if ($automake && !@$dnafiles && $outformats{'fasta'}) {
$self->preMake('dna');
$dnafiles = $self->getFiles('dna', $chromosomes);
}
# if ($automake && !@$featfiles && ($outformats{'fasta'} || $outformats{'gnomap'})) { #(grep /fasta|gnomap/, @outformats))
# $self->preMake('fff'); # not active?
# $featfiles = $self->getFiles('fff', $chromosomes);
# }
if ($DEBUG) {
my @cn= @$chromosomes; # FIXME too long to list for some, no_csomesplit
@cn= (@cn[0..20],".. more ..") if(@cn>30);
print STDERR "make chromosomes= @cn\n";
my @fn= map {$_->{name}} @$featfiles; print STDERR "with featfiles= @fn\n";
my @dn= map {$_->{name}} @$dnafiles; print STDERR "with dnafiles= @dn\n";
}
if (delete $outformats{'fasta'}) {
$self->missingData( $featfiles, "fff", "fff files; make with -format fff");
$self->missingData( $dnafiles, "dna", "dna files; make with -dnadump");
my $writer= $self->getWriter('fasta');
push @results, $writer->makeFiles(%args,
infiles => $featfiles, chromosomes => $chromosomes);
}
if (delete $outformats{'blast'}) {
my $fafiles = $self->getFiles( 'fasta', $chromosomes);
$self->missingData( $fafiles, 'fasta', "fasta files; make with -format fasta");
my $writer = $self->getWriter('blast'); # this works; eval new writer
push @results, $writer->makeFiles( %args,
infiles => $fafiles, chromosomes => $chromosomes );
}
if (delete $outformats{'gnomap'}) { # move to generic handler below
$self->missingData( $featfiles, "fff","fff files; make with -format fff");
my $writer= $self->getWriter('gnomap');
push @results, $writer->makeFiles(%args,
infiles => [ @$featfiles, @$dnafiles ], chromosomes => $chromosomes);
# needs $featfiles
}
# my @moreformats= grep !/(fff|gff|dna|fasta|blast|gnomap)/,@outformats;
my @moreformats= grep { $outformats{$_} } @outformats; # preserve call-order
foreach my $fmt (@moreformats) {
my $writer= $self->getWriter($fmt);
my $fset = $self->getFilesetInfo($fmt);
my $infiles;
if($fset->{input} =~ /feature_table/) {
$infiles = $self->getFiles('feature_table', $chromosomes);
$self->missingData( $infiles, 'feature_table', "feature_table files; make with -featdump");
} else {
# ($fset->{input} =~ /fff/)
$infiles= [ @$featfiles, @$dnafiles ]; #?? always
# $self->missingData( $featfiles, "fff", "fff files; make with -format fff");
# $self->missingData( $dnafiles, "dna", "dna files; make with -dnadump");
}
if (!$writer) { warn "no writer for $fmt\n"; }
else { push @results, $writer->makeFiles( %args,
name => $fmt,
filesetinfo => $fset,
infiles => $infiles,
chromosomes => $chromosomes);
}
}
=item more formats using SQL queries
-- need a Writer to handle dumpFeatures : see BulkWriter now
if (delete $outformats{'go_association'}) {
my $fset = $self->getFilesetInfo('go_association');
if($fset) {
# my ($ovfiles,$ovmissing)= $self->getDumpFiles(['go_association'], $fset);
my $ovfiles = $self->dumpFeatures($fset, undef, "colnames");
my $ovlist= join(" ",map {$_->{name}} @$ovfiles);
push @results, "output:$ovlist";
}
}
=cut
$self->gzipFiles( \@outformats, $chromosomes );
my $lok= $self->makeCurrentLink()
if (@results && $self->config->{make_current});
#? put most of above in eval{} block so we return error info if failed ??
push(@results,"\nBulkfiles are located at ".$self->getReleaseDir())
if(@results);
return join(", ",@results); #what?
}
=item writeDocs( $docs or $self->config->{doc})
print docs from config file
.. move this into own BulkWriter subclass ?
need some fix to writeDocs for doc->path at top level or not-releasedir
=cut
sub writeDocs
{
my ($self, $docs)= @_;
my $ndoc= 0;
if (ref $docs) {
# check for 1 or many (name keys, darn xmlsimple); is tag 'name' or 'id' ?
my $reldir = $self->getReleaseDir();
my $datadir= $self->config->{datadir}; # must exist
my $species= $self->speciesFull();
if ($docs->{content}) {
my %dd= ();
if ($docs->{name}) { %dd= ( $docs->{name} => $docs ); }
elsif ($docs->{id}) { %dd= ( $docs->{id} => $docs ); }
else { %dd= ( 'untitled' => $docs ); }
$docs= \%dd;
}
foreach my $dname (sort keys %$docs) {
next if ($docs->{$dname}->{hidden}); ## == 1
my $data = $docs->{$dname}->{content} || '';
my $dpath= $docs->{$dname}->{path} || $dname;
$dpath =~ s/\${datadir}/$datadir/;
$dpath =~ s/\${species}/$species/;
my $norel=($dpath =~ m/$datadir/ || $dpath =~ m,^/,);
my $fn= ($norel) ? $dpath : catfile( $reldir, $dpath);
print STDERR "write doc $dname $fn\n" if $DEBUG;
if( open(DOC,">$fn") ) { print DOC $data; close(DOC); $ndoc++; }
else { warn "ERROR: cant write $fn\n"; }
}
}
print STDERR "writeDocs n=$ndoc\n" if $DEBUG; # debug
return $ndoc;
}
=item getWriter
Replaced getXxxWriter with generic getWriter('type')
=cut
sub getWriter
{
my ($self, $type)= @_;
my $finfo= $self->getFilesetInfo($type);
if (ref $finfo && $finfo->{handler}) {
my $configfile= $finfo->{config};
my $pkg= $finfo->{handler};
if($pkg eq "Bulkfiles") { return $self; } #?? can we do this; needs BulkWriter methods
unless($pkg =~ /\:\:/) { $pkg= "Bio::GMOD::Bulkfiles::".$pkg; }
my $eval=
"use $pkg;
$pkg->new( configfile => \$configfile, fileinfo => \$finfo,
handler => \$self, debug => \$DEBUG, showconfig => \$self->{showconfig},
);";
##print STDERR "getWriter: eval $eval\n" if $DEBUG;
my $writer= eval $eval;
if ($@) { warn $@; die if($self->{failonerror}); }
return $writer if ref $writer;
}
# print STDERR "getWriter('$type'): eval failed\n" if $DEBUG;
# ## old way
# CASE: {
# $type eq 'fasta' && return $self->getFastaWriter();
# $type eq 'blast' && return $self->getBlastWriter();
# $type eq 'fff' && return $self->getFeatureWriter();
# $type eq 'gff' && return $self->getFeatureWriter();
# $type eq 'gnomap' && return $self->getGnomapWriter();
# }
warn "no writer module for $type\n";
}
#===================================================
=item getChromosomeTable
locate and read feature dump of chromosomes (or equivalent parts)
2L 1 22217931 0 chromosome_arm 2L 2L 1 species Drosophila_melanogaster
2R 1 20302755 0 chromosome_arm 2R 2R 2 species Drosophila_melanogaster
3L 1 23352213 0 chromosome_arm 3L 3L 4 species Drosophila_melanogaster
3R 1 27890790 0 chromosome_arm 3R 3R 3 species Drosophila_melanogaster
4 1 1237870 0 chromosome_arm 4 4 5 species Drosophila_melanogaster
U 1 11561901 0 chromosome_arm U U 7 species Drosophila_melanogaster
X 1 21780003 0 chromosome_arm X X 6 species Drosophila_melanogaster
oct04 -- add option to create 'chr-U' from bag of Unknown_* golden_path entries
in getChromosomeTable. Have some 2000+ Unknown_group and Unknown_singleton
golden_path/ultra_scaffold entries in Dpse r2.1.
need also map contig -> contig_contig -> ultra_scaffold
Contig1083_Contig6433 0 11665 1 3 contig Contig1083 3209332 parent_oid 3208389:1
Unknown_group3 0 15476 1 3 golden_path_region Contig1083_Contig6433 Contig1083_Contig6433 3208389 parent_oid 3798908:1
return hash {
'3R' => {
'length' => 27890789,
'start' => 1,
'oid' => 3,
'strand' => '0',
'type' => 'chromosome_arm',
'name' => '3R'
},
4 => {
'length' => 1237869,
'start' => 1,
'oid' => 5,
'strand' => '0',
'type' => 'chromosome_arm',
'name' => 4
},
}
=cut
sub getChromosomeTable
{
my $self= shift;
my $config= $self->config;
if (defined $config->{chromosome}) { return $config->{chromosome}; }
my $chromosome= {};
my $chrparts = {}; # for dpse map Unknown ultra_scaffold/golden_path_fragment to U
my $part2chr = {}; # map golden_path_region name to chr-arm name
my $chrpartpattern= $config->{chrpart_pattern};
my %orgset = ();
my $nozombiechromosomes= $config->{nozombiechromosomes};
# dpse chado duplicate 0-length chromosome entries
## need remapArm() here, but csomeset is result of this method ...
## BUT featdump tables are not remapped until csomeSplit(); cant use remapped names til after that
#? allow only one species per make run ?
my $myspp= $config->{species}; $myspp =~ s/ /_/g;
my $myorg= $config->{org} || $self->speciesAbbrev($myspp);
my $fileset= $self->getDumpFiles(['chromosomes']);
my $path= (ref $fileset) ? $fileset->[0]->{path} : undef;
if ($path && open(CF,$path)) {
while() {
next unless(/^\w/);
next if(/^arm\tfmin/); # header from sql out -- should be 'chromosome' or 'chr' instead of 'arm'
chomp;
my $oldarm;
my ($arm,$fmin,$fmax,$strand,$orgid,$type,$name,$id,$oid,$attr_type,$attribute)
= split("\t");
next unless($id); #?
## sgdlite uses messy chr ID -- use name instead here ? better: $arm is best of both
## fb-dpse has zero-length chromosomes, and duplicates (some zero some not)
## ? skip zero len csomes ? for output at least
## use $strand == rank here -- assume input file is ordered by that.
## need ($arm,$golden_path,...)= mapChr($arm)
## ? need some compound chr{arm} with multiple ids for Unknown bag?
($arm,$fmin,$fmax,$strand,$oldarm)= $self->remapArm($arm,$fmin,$fmax,$strand,undef);
next unless($arm && $arm ne 'skip'); #?
my $species= ($attr_type eq 'species') ? $attribute : $config->{species};
$species =~ s/ /_/g;
my $org= $self->speciesAbbrev($species);
$org= $self->speciesAbbrev($orgid) unless($org);
$org= "null" unless($org);
next unless($myorg eq $org || $myspp eq $species); #?? dec05; for dpse+dmel db
next if ($nozombiechromosomes && $fmax <= $fmin);
my $org4= $self->speciesAbbrev($orgid, "org4");
my $chrvals= {
arm => $arm,
oldarm => $oldarm, # rarely differs from arm ...
name => $name || $id,
id => $id,
type => $type,
start => $fmin,
length => ($fmax - $fmin + 1),
strand => $strand,
oid => $oid,
species => $species,
orgid => $orgid, ## NEED THIS, now all feature_table have organism_id
org => $org,
};
my($genus,$spp)= split(/_/,$species,2);
$orgset{$orgid}= {
organism_id => $orgid,
abbreviation => lc($org),
org4 => $org4,
genus => $genus,
species => $spp,
fullspecies => $species,
};
if ($chrpartpattern && $type =~ /$chrpartpattern/) { #== golden_path_fragment|golden_path_region, other?
unless($chrparts->{$arm}) { $chrparts->{$arm}= []; }
push(@{$chrparts->{$arm}}, $chrvals);
$part2chr->{$id}= $arm;
}
else {
#? check type =~ ^(chromosome_arm|golden_path|ultra_scaffold)$
$chromosome->{$arm}= $chrvals;
}
}
close(CF);
}
$chromosome->{'_part2chr'}= $part2chr if( %$part2chr );
foreach my $arm (keys %$chrparts) {
unless($chromosome->{$arm}) { ## make pseudochr from parts?
$chromosome->{$arm}= {
arm => $arm, name => $arm, id => $arm,
type => 'golden_path', # fixme
start => 1,
length => 0, # fixme
strand => 0, oid => 0, species => '', org => '', # fixme
};
}
$chromosome->{$arm}->{parts}= $chrparts->{$arm};
}
my $organisms = $self->getOrganismTable(\%orgset);
$config->{chromosome}= $chromosome;
$self->getChromosomes();
warn "N chromosomes=",scalar(keys %{$chromosome}),"\n" if $DEBUG;
return $chromosome;
}
sub getOrganismTable
{
my $self= shift;
my($orgset)= @_;
## add read from overview tables/organisms-overview.txt
## flds: qw(Organism_id Abbreviation Genus Species Common_name N_features Comment)
unless($self->{did_organism}) {
my $overviewset = $self->getFilesetInfo('overview');
my $tabdir = $self->getReleaseSubdir( $overviewset->{path} || 'tables/');
my $tabfile= catfile( $tabdir, "organisms-overview.txt");
my @colheads;
if( open(DOC,$tabfile) ) {
$orgset= {} unless(ref $orgset);
$self->{did_organism}++;
while(){
chomp; my @col= split "\t";
if(/^Organism_id/i){ @colheads= @col; }
elsif(/^\d/ && scalar(@col) > 3) {
foreach (@col) { $_='' if($_ eq "\\N"); }
my($orgid,$abbreviation,$genus,$species,$common,@xtra)= @col;
my $fullspecies= $genus."_".$species;
unless($abbreviation) { $abbreviation= $self->speciesAbbrev($fullspecies); }
my $org4= $self->species4letter($fullspecies);
$orgset->{$orgid}= {
organism_id => $orgid,
org4 => $org4,
abbreviation => lc($abbreviation),
genus => $genus,
species => $species,
fullspecies => $fullspecies,
from_db => 1,
};
}
}
close(DOC);
}
}
my $organisms= $self->config->{organism};
$organisms= {} unless(ref $organisms);
my $norgs=0;
if(ref $orgset) {
foreach my $orgid (keys %$orgset) {
my $orgref= $orgset->{$orgid};
$organisms->{$orgid}={} unless(ref $organisms->{$orgid});
next if( $organisms->{$orgid}->{from_db} && ! $orgref->{from_db});
$organisms->{$orgid}->{organism_id}= $orgid;
my $fullspecies= $orgref->{fullspecies} || $orgref->{genus}."_".$orgref->{species};
$organisms->{$orgid}->{fullspecies}= $fullspecies;
$organisms->{$orgid}->{species}= $fullspecies; #?? need this same as fullspecies??
$organisms->{$orgid}->{from_db}= $orgref->{from_db} || 0;
$organisms->{$orgid}->{org4}= $self->species4letter($fullspecies);
my $abbrev= $orgref->{abbreviation} || $self->speciesAbbrev( $fullspecies );
$organisms->{$orgid}->{abbreviation}= $abbrev;
$organisms->{$abbrev}= $organisms->{$orgid}; # copy for orgid lookup
$norgs++;
}
warn "Organisms n_entries=$norgs\n" if $DEBUG;
}
$self->config->{organism}= $organisms;
return $organisms;
}
sub getChromosomes
{
my $self= shift;
my $config= $self->config;
unless(ref $config->{chromosomes}) {
my $chromosome= $self->getChromosomeTable();
my @csomes= grep !/^_/, sort keys %$chromosome;
$config->{chromosomes}= \@csomes;
}
return $config->{chromosomes};
}
# this needs to be a configuration choice (how many genus,species letters)
# and combine with chado table abbreviation, locase or not, ... have too many
# abbreviations now, used in file names. Need simpler common abbrev.
sub species4letter
{
my $self= shift;
my ($spp)= @_;
$spp ||= $self->config->{species} || $self->config->{org};
$spp =~ s/ /_/g;
my $spattern= $self->config->{species_short_pattern}
|| '^(\w)[^_]*_(\w{1,3})';
my( $ga, $sa) = $spp =~ m/$spattern/;
return lc($ga.$sa) if($ga and $sa); #? should lc() be a _pattern option?
return lc( substr($spp,0,4) ); #??
}
sub speciesAbbrev
{
my $self= shift;
my ($spp, $org4letter)= @_;
$spp = $spp || $self->config->{species} || $self->config->{org} || "";
$spp =~ s/ /_/g;
my $organisms= $self->config->{organism};
if (ref $organisms) {
local $^W=0; # kill warnings of undef orgset values
foreach my $org (reverse sort keys %{$organisms}) {
my $orgset= $organisms->{$org};
if ($spp eq $org
|| $spp eq $orgset->{fullspecies}
|| $spp eq $orgset->{species}
|| $spp eq $orgset->{organism_id}
) {
if($org4letter) {
my $org4= $orgset->{org4} || $self->species4letter($spp);
return $org4;
}
elsif ($org =~ /\d+/) { # watchout for org == orgid here
my $abbr= $orgset->{abbreviation};
return $abbr if $abbr;
}
else { return $org; }
}
}
}
return $self->species4letter($spp);
}
## fix this; should always prefer $self->config->{species}; ?
## $organisms->{$orgid}->{fullspecies} second choice?
## problems with ' ',_ in species such as
## genus=Anopheles species='gambiae str. PEST'
sub speciesFull
{
my $self= shift;
my ($org)= @_;
my $species= '';
$species= $self->config->{species};
$species =~ s/ /_/g;
return $species if($species =~ m/_/); # problems
$org= $org || $species || $self->config->{org};
$org= $self->speciesAbbrev($org);
my $organisms= $self->config->{organism};
if(ref $organisms && $organisms->{$org}->{species}) {
$species= $organisms->{$org}->{species};
$self->config->{species}= $species if ($species =~ m/_/);
}
return $species;
}
sub speciesFull_OLD
{
my $self= shift;
my ($org)= @_;
my $species= '';
unless($org) {
$species= $self->config->{species};
$species =~ s/ /_/g;
return $species if($species =~ m/_/);
}
$org= $org || $species || $self->config->{org};
$org= $self->speciesAbbrev($org);
my $organisms= $self->config->{organism};
$species= $organisms->{$org}->{species} if ($organisms->{$org}->{species});
$self->config->{species}= $species if ($species =~ m/_/);
return $species;
}
=item splitFFF
split flat feature format line to fields
($type,$name,$cytomap,$baseloc,$id,$dbxref,$notes,$chr)
= splitFFF($fffeature, $chr)
=cut
sub splitFFF
{
my( $self, $fffeature, $chr)= @_;
my($type,$name,$cytomap,$baseloc,$id,$dbxref,$notes,$bstart);
chomp($fffeature);
my @v= split "\t", $fffeature;
foreach (@v) { $_='' if $_ eq '-'; }
my $ffformat = $self->{ffformat} || 0; #? test always
unless( $ffformat > 0 ) {
if ( @v > 7 || ($v[0] =~ /^\w/ && $v[1] =~ /^\d+$/)) { $ffformat= 2; }
else { $ffformat= 1; }
}
if ($ffformat == 1) { ($type,$name,$cytomap,$baseloc,$id,$dbxref,$notes)= @v; }
elsif ($ffformat == 2) { ($chr,$bstart,$type,$name,$cytomap,$baseloc,$id,$dbxref,$notes)= @v; }
$self->{gotffformat}= $ffformat;
return ($type,$name,$cytomap,$baseloc,$id,$dbxref,$notes,$chr);
}
=item $interfff= intergeneFromFFF2($chr,$fff1,$fff2)
return fff line for region between two features
(any feature types ok, intergenic regions most interesting)
=cut
sub intergeneFromFFF2
{
my ($self, $chr, $fff1, $fff2)= @_;
my $newfff='';
my ($type,$name1,$cytomap,$baseloc1,$id1,$dbxref,$notes);
my ($chr2,$name2,$baseloc2,$id2);
my $intergenetype='intergene'; # FIXME
($type,$name1,$cytomap,$baseloc1,$id1,$dbxref,$notes,$chr)
= $self->splitFFF( $fff1, $chr);
my($start1,$stop1,$strand1)= $self->maxrange($baseloc1);
($type,$name2,$cytomap,$baseloc2,$id2,$dbxref,$notes, $chr2)
= $self->splitFFF( $fff2, $chr);
my($start2,$stop2,$strand2)= $self->maxrange($baseloc2);
if ($chr eq $chr2 && ($stop1 + 2) < $start2) {
my $iname= "$name1/$name2";
my $iid= "$id1/$id2";
my $interloc= ($stop1+1)."..".($start2-1);
if ($self->{gotffformat} == 2) {
$newfff= join("\t", $chr, $stop1, $intergenetype, $iname,'-', $interloc, $iid,$dbxref,$notes);
} else {
$newfff= join("\t", $intergenetype, $iname,'-', $interloc, $iid,$dbxref,$notes);
}
}
return $newfff;
}
## no_csomesplit: replace dna/ files w/ chado db calls for all?? see dumpChromosomeBases
sub getDnaSeqFromDb
{
my($self,$chr)= @_;
my $csomeset= $self->getChromosomeTable();
my $id= $csomeset->{$chr}->{id} || $chr;
my $bases= $self->getBasesFromDb($id,1);
unless ($bases) {
warn "getDnaSeqFromDb: no bases for segment:$chr\n";
die if $self->{failonerror};
}
return $bases;
}
sub getDnaSeq
{
my ($self, $chr)= @_;
my $seq= $dnaseqs{$chr};
unless($seq) {
my $no_csomesplit= $self->config->{no_csomesplit} || 0; # 2007oct: old default, change for many small scaffolds
if($no_csomesplit) {
my $dna= $self->getDnaSeqFromDb($chr);
$seq= Bio::PrimarySeq->new( -id=>$chr, -seq => $dna);
print STDERR "getDnaSeqFromDb $chr, length=",$seq->length(),"\n" if $DEBUG;
} else {
my $dnafile= $self->dnafile($chr); #"$dnadir/dna-$chr.raw";
$seq= Bio::GMOD::Bulkfiles::MyLargePrimarySeq->new( -id=>$chr, -file => $dnafile);
print STDERR "open dnafile $dnafile, length=",$seq->length(),"\n" if $DEBUG;
}
$dnaseqs{$chr}= $seq;
}
return $seq;
}
sub getBases
{
my($self, $usedb,$type,$chr,$baseloc,$id,$name,$subrange,$makeaa)= @_;
my $bases= undef;
if($usedb && $id) {
$bases= $self->getBasesFromDb($id);
use constant FIX_AAMISSING_STOP => 1;
if(FIX_AAMISSING_STOP) {
if($makeaa && $bases) {
$bases .='*' unless($bases =~ /\*$/); # annoying non-standard store
if(ref $self->{gotbases}) { $self->{gotbases}->{'md5checksum'}=0; } # no longer valid
}
}
return $bases if($bases || $self->{failonerror});
# >> no; try from files: $self->{skiponerror});
}
unless ($bases) {
$bases= $self->getBasesFromFiles($type,$chr,$baseloc,$name,$subrange);
if($makeaa) { $bases= $self->dna2protein($bases); }
}
return $bases;
}
# but see Bio/GMOD/DB/Config.pm
sub dbiDSN
{
my ($self, $dsn)= @_;
my $config= $self->config;
my ($dbuser,$dbpass)=("","");
if ($dsn && $dsn =~ /^dbi:/) { $self->{dsn}= $dsn; }
if (ref $config->{db}) {
my $dbname= $config->{db}->{name};
my $relid= $config->{relid}; #? or now use promoted {release_dbname}
my $reldb= ($relid && defined $config->{release}->{$relid})
? $config->{release}->{$relid}->{dbname} :'';
$dbname= $reldb if ($reldb);
unless($dbname) { warn "missing dbname"; die if $self->{failonerror}; }
## ? handle dbi:mysql:database=dmel_r41_20050207;host=localhost;port=3306;mysql_socket=/tmp/fbmysql.sock
if($config->{db}->{dsn}){
$dsn = $config->{db}->{dsn};
} else {
$dsn = "dbi:" . $config->{db}->{driver} || "Pg";
$dsn .= ":dbname=" .$dbname;
$dsn .= ";host=" .$config->{db}->{host} if $config->{db}->{host};
$dsn .= ";port=" .$config->{db}->{port} if $config->{db}->{port};
#?? $self->{dsn}= $dsn;
}
$dbuser= $config->{db}->{user} if $config->{db}->{user};
$dbpass= $config->{db}->{password} if $config->{db}->{password};
}
## if ($self->{dsn}) { $dsn= $self->{dsn}; }
return (wantarray) ? ($dsn,$dbuser,$dbpass) : $dsn;
}
sub dbiConnect
{
my ($self)= @_;
my $dbh= $self->{dbh};
unless($dbh) {
my $tdsn= $self->dbiDSN();
print STDERR "DBI->connect( $tdsn )\n" if $DEBUG;
$dbh = DBI->connect( $self->dbiDSN() )
or die("unable to open db( $tdsn )"); # throw ?
$self->{dbh}= $dbh;
}
return $dbh;
}
=item getSeqSql($sqlconf)
read in config file with feature dump sql scripts
=cut
sub getSeqSql
{
my ($self, $sqlconf, $sqlenv)= @_;
$sqlconf = 'chadofeatsql' unless($sqlconf);
$sqlenv= $self->config unless (ref $sqlenv);
# print STDERR "sqlenv.$sqlconf: ",join("\n ", map{ $_."=".$sqlenv->{$_}} keys %$sqlenv ),"\n"
# if ($DEBUG>1);
my $config2= $self->{config2}; #?? Config2 object, not hash
my $seqsql = $self->{$sqlconf} || '';
unless($seqsql) {
$seqsql= $config2->readConfig( $sqlconf, {Variables => $sqlenv}, {} );
print STDERR $config2->showConfig($seqsql, { debug => $DEBUG })
if ($self->{showconfig} && $DEBUG>1);
}
# has undefined Variables? try to define
my $sqxml= $config2->showConfig( $seqsql, { debug => 0 });
if ( $sqxml =~ m/\$\{/ ) {
eval { $seqsql= $self->{config2}->updateVariables( $seqsql, { Variables => \%ENV, replace =>1 }, ); };
if ($@) { warn "Config2: getSeqSql=$sqlconf; err: $@"; }
}
if ( $sqxml =~ m/\$\{/ && ref($seqsql->{ENV_default}) ) {
my %env= %{$seqsql->{ENV_default}};
foreach my $k (keys %env) { $env{$k}= $sqlenv->{$k} if($sqlenv->{$k}); }
$seqsql= $config2->readConfig( $sqlconf, {Variables => \%env}, {} );
}
# print STDERR "seqsql.$sqlconf xml: ", $config2->showConfig( $seqsql, { debug => 0 }),"\n"
# if ($DEBUG>1);
$self->{$sqlconf}= $seqsql;
return $seqsql;
}
=item validateVariables
validate config values for seq_ontology, golden_path, species??
to exist in database
=cut
sub validateVariables
{
my ($self)= @_;
return if $self->{didvalidatevars};
$self->{didvalidatevars}= 1;
warn "Validating...\n" if($DEBUG||$self->{automake});
my $note="";
my $doinspect= ($DEBUG or $self->{verbose} or $self->config->{newuser}) ? 1 : 0;
## ? newuser, verbose, DEBUG, ...
my $fdump = $self->config->{featdump};
my $seqsql = $self->getSeqSql($fdump->{config},$fdump->{ENV});
my $sqltag = $fdump->{tag} || "feature_sql";
my $sqltype = 'validate'; # other types?
my @sqlparam=();
## tag names here are now seq_ontology_check, golden_path_check
my $dbh= $self->dbiConnect();
my @targets= sort keys %{$seqsql->{$sqltag}};
foreach my $tgname (@targets) {
my $fs = $seqsql->{$sqltag}->{$tgname};
my $type= $fs->{type};
my $sql = $fs->{sql};
next unless($sql && ( $type =~ m/\b$sqltype\b/) );
$note.= "$tgname = ";
print STDERR "$type:$tgname\n" if $DEBUG; # sql=$sql
# local $^W=0; # kill warnings of undef values
my $nrow= 0;
my $err="";
my $sth = $dbh->prepare($sql) or $err= "Failed to prepare $sql";
unless($err) { $sth->execute() or $err= "Failed to execute sql"; }
unless($err) { while (my $row = $sth->fetchrow_arrayref) { $nrow++; } }
$sth->finish();
my $sql_inspect = $fs->{sql_inspect}; # optionally show this inspect even if no err
if($err or $nrow == 0) {
$err .= "*** Failed validation check for $tgname\n";
$err .= $fs->{warning}."\n" if($fs->{warning});
$err .= "$tgname sql=$sql\n";
warn $err; $err= "Invalid $tgname\n";
$doinspect= 1;
}
if($doinspect and $sql_inspect) {
##??? problems with this for chromosome checks (using a join) very long for parameciumDB
## why? is it bug in sql/postgres/database tuning? using some scan instead of index search?
## explain high cost in chado pg is with Sort for Group
## Sort (cost=178872.73..179239.36 rows=146652 width=560); Sort Key: f.type_id, t.name
warn "\nInspecting database for $tgname values...\n";
print STDERR ("-") x 60, "\n";
my $result= $self->getFeaturesFromDb( *STDERR, $sql_inspect, undef, "colnames");
print STDERR ("-") x 60, "\n";
}
if ($err) {
warn $err; die if($self->{failonerror}); #?? ||$self->{automake}
$note.="Invalid ";
}
else { $note.="Ok "; }
}
warn "Variable checks: $note\n" if($DEBUG||$self->{automake});
}
=item updateSqlViews
add views to db used by sql feature calls.
this may fail if one lacks update permissions; assume
user knows about such
=cut
sub updateSqlViews
{
my ($self, $seqsql, $sqltag)= @_;
return if $self->{didsqlviews};
$self->{didsqlviews}= 1;
unless($seqsql) {
my $fdump = $self->config->{featdump};
$seqsql = $self->getSeqSql($fdump->{config},$fdump->{ENV});
}
$sqltag ||= "feature_sql";
my $sqltype = 'view'; # other types? procedures ?
my $dbh= $self->dbiConnect();
my @targets= sort keys %{$seqsql->{$sqltag}};
foreach my $tgname (@targets)
{
my $fs= $seqsql->{$sqltag}->{$tgname};
my $type= $fs->{type};
my $sql = $fs->{sql};
unless($sql && ( $type =~ m/\b$sqltype\b/) ) { next; }
print STDERR "do sql $tgname $type\n" if $DEBUG;
my $result = $dbh->do($sql) or warn "unable to do sql $tgname $type";
}
}
=item dumpFeatures
dumpFeatures($fdump, $sqlconf) - extract feature_table s from chado db using
feature sql config info
-- add other config items for sql dumps - organism_table; lists; ..
-- use fileset instead of featdump
=cut
sub dumpFeatures
{
my ($self, $fdump, $sqlconf, $dumpflags)= @_;
my @files=();
$dumpflags ||=""; ## add colnames for overview.txt; look in $fdump->{dumpflags} ??
$fdump ||= $self->config->{featdump};
$sqlconf ||= $fdump->{config_sql} || $fdump->{config};
my $seqsql = $self->getSeqSql($sqlconf,$fdump->{ENV});
$self->updateSqlViews($seqsql, $fdump->{tag});
my $outpath= $self->getReleaseSubdir( $fdump->{'path'} || catdir( $self->config->{TMP}, "featdump") );
$fdump->{'path'}= $outpath; # save for reuse
my $sqltag = $fdump->{tag} || "feature_sql";
my $sqltype = $fdump->{type};
my $targets = $fdump->{target}; # should be array ?
unless($targets) { my @tg= sort keys %{$seqsql->{$sqltag}}; $targets= \@tg; }
unless(ref $targets) { $targets= [ $targets ]; }
foreach my $tgname (@$targets)
{
my $fs= $seqsql->{$sqltag}->{$tgname};
unless($fs) { warn "no sql dump target $tgname in $sqlconf"; next; }
my $type= $fs->{type};
my $sql = $fs->{sql};
my $outn= $fs->{output} || $tgname.".tsv";
unless($sql && (!$sqltype || $type =~ m/\b$sqltype\b/)) {
warn "dumpFeatures skip: $tgname/$type<>$sqltype\n";
next; } #??
my $outf= catfile($outpath,$outn);
my $outh= new FileHandle(">$outf");
print STDERR "sql dump $tgname $type $outf\n" if $DEBUG;
my $nout= $self->getFeaturesFromDb( $outh, $sql, undef, $dumpflags);# \@sqlparam ?
print STDERR "sql dump $tgname n rows=$nout\n" if $DEBUG;
close($outh);
my $fixme = $fs->{script}; # may be array/hash of scripts ?
if ($fixme && $fixme->{type} eq 'postprocess') {
my $shell= $fixme->{shell} || $fixme->{language};
my $spath= catfile( $self->config->{TMP}, $fixme->{name});
my $fixinput= $outf;
print STDERR "postprocess $shell $spath $fixinput\n" if $DEBUG;
open(SH,">$spath"); print SH $fixme->{content}; close(SH);
my $sresult= `$shell $spath $fixinput`; #?? what of perl params
# how do i pipe table into script ?? and out again to replace old
# this works: perl -i.old rdump $r/tmp/featdump/analysis.tsv
}
## changed keys: name=> target; file=> name
push(@files, { path => $outf, type => $type, target => $tgname, name => $outn, });
}
return \@files;
}
=item getFeaturesFromDb
$nrows = self->getFeaturesFromDb( $outh, $sql, \@sqlparam, $flags)
$table = self->getFeaturesFromDb( undef, $sql, \@sqlparam, $flags)
=cut
sub getFeaturesFromDb
{
my ($self, $outh, $sql, $sqlparam, $flags)= @_;
my $err="";
my $dbh = $self->dbiConnect();
my $sth = $dbh->prepare($sql) or $err= "failed to prepare $sql";
if (ref $sqlparam) { $sth->execute(@$sqlparam) or $err= "failed to execute sql" ; }
else { $sth->execute() or $err= "failed to execute sql"; }
if ($err) { warn $err; die if($self->{failonerror}); return -1; }
## analysis sql gets lots of: Use of uninitialized value in join
## OUT OF MEM error here with ParameciumDB during fetchrow to tmp/featdump file
## get same memerr when run psql -f features.sql, so it isn't Bulkfiles problem !
## *** malloc: vm_allocate(size=67108864) failed (error code=3); could be gffattr_gmodel
local $^W=0; # kill warnings of undef values
my $result= 0;
if($outh) {
print $outh join("\t",@{$sth->{NAME}})."\n" if($flags=~/colname/i);
while (my $row = $sth->fetchrow_arrayref) {
print $outh join("\t",@$row),"\n"; $result++;
}
} else {
$result = join("\t",@{$sth->{NAME}})."\n" if($flags=~/colname/i); # is this ok?
while (my $row = $sth->fetchrow_arrayref) {
$result .= join("\t",@$row)."\n";
}
}
$sth->finish();
return $result;
}
=item $bases= getBasesFromDb( $uniquename)
=cut
sub getBasesFromDb
{
my ($self, $uniquename, $iscsome)= @_;
my $dbh= $self->dbiConnect();
my $sql="";
$sql= $self->config->{dnadump}->{sql_csome} if ($iscsome);
$sql ||= $self->config->{dnadump}->{sql};
$sql ||= "select feature_id, residues, md5checksum, seqlen, name from feature where uniquename = ?";
my $err="";
my $sth = $dbh->prepare($sql) or $err="unable to prepare feature_id";
unless($err) { $sth->execute($uniquename) or $err="failed to execute feature_id"; }
if ($err) { warn $err; die if ($self->{failonerror}); return undef; }
## check for >1 row! e.g. flybase dpse has 2+ csome entries, 1 w/o bases, for same ID
my($hashref,$feature_id,$bases)=(undef,undef,"");
while (my $nextrow = $sth->fetchrow_hashref) {
$hashref = $nextrow;
$feature_id = $$hashref{'feature_id'};
$bases= $$hashref{'residues'} || "";
last if $bases;
}
$self->{gotbases}= $hashref; ## return md5checksum, seqlen fields ;stick all other fields into self
$sth->finish();
print STDERR "getBasesFromDb $uniquename -> $feature_id ; n bases=",length($bases),"\n"
if (!$bases || $DEBUG > 1);
return $bases;
}
sub getLastBasesFeature
{
my ($self)= @_;
return $self->{gotbases};
}
sub maxrange
{
my ($self, $range)= @_;
my ($pre, $suf,$start,$stop, $b, $u);
$start= -9; $stop= $start;
$range =~ s/^([^\d<>-]*)//; $pre= $1;
$range =~ s/(\D*)$//; $suf= $1;
my $strand= ($pre =~ /^complement/) ? -1 : 1;
if ($range =~ m/^([<>]*)([\d\-]+)/) { $u= $1; $start= $2; $start-- if ($u eq '<'); }
if ($range =~ m/([<>]*)([\d\-]+)$/) { $u= $1; $stop= $2; $stop++ if ($u eq '>'); }
return ($start,$stop,$strand);
}
sub dna2protein
{
my($self, $dna)= @_;
my $aa= undef;
eval {
require Bio::Tools::CodonTable;
my $codonTable = Bio::Tools::CodonTable->new();
$aa = $codonTable->translate($dna);
};
if ($@) {
warn "dna2protein: err: $@";
die if ($self->{failonerror} || $self->{automake});
}
return $aa;
}
sub getBasesFromFiles
{
my($self, $type,$chr,$baseloc,$name,$subrange)= @_;
my $bases= undef;
my $gotloc= $baseloc;
my $dnaseq= $self->getDnaSeq($chr);
if ($dnaseq) {
my ($start,$stop,$strand)= $self->maxrange($baseloc);
my ($subrb,$subre,$rs)=(0,0,0);
if ($subrange) {
# need some more logic in $subrange to get just upstream or downstream sections
# readseq's: start +/- offset1, stop +/- offset2; eg. (start-2000,start); (stop,stop+2000); (start-2000,stop); (start,stop+2000)
my $maxseq= $dnaseq->length();
($subrb,$subre,$rs)= $self->maxrange($subrange);
if ($subrb) { $start += $subrb ; $start=1 if $start<=0; } # need dnaseq min/max !
if ($subre) { $stop += $subre ; $stop=$maxseq if $stop> $maxseq; } # need dnaseq min/max !
}
my $range= $baseloc;
$range =~ s/^[\w]*\(?//; ##s/^([^\d<>-]*)//;
$range =~ s/\)?\s*$//;
my @locs= split(/,/,$range);
=item handle Bioperl version -strand effects
080620: Bioperl versions are unpredictable in how they handle
multiexon -strand thru Location::Split and (Large)PrimarySeq use of revcomp for -split locs.
Avoid this by feeding in only +strand exons, but for rare cases of mixed strand genes,
then do explicit revcom on final CDS-dna when needed.
Input baseloc is simple form: complement(x,y,z), NOT messy join(comp(x),comp(y),comp(z)).
add_sub_Location(start,stop,exon+strand) will always use exon+strand even when CDS is -topstrand.
This still allows per exon -strand, e.g., compl(x,y,compl(z)) => (-x, -y, +z)
=cut
if (@locs>1) {
## avoid bad strand() handling in Bio::Location::Split();
my $sloc = Bio::GMOD::Bulkfiles::MySplitLocation->new();
my $topstrand= $strand;
# $sloc->strand($topstrand); # Bioperl action here was bad; AVOID THIS
# if ($topstrand < 0) { @locs = reverse @locs; } # not right now; AVOID THIS
foreach my $loc (@locs) {
($start,$stop,$strand)= $self->maxrange($loc);
# $strand= -$strand if ($topstrand < 0); # AVOID THIS
$sloc->add_sub_Location( new Bio::Location::Simple(
-start => $start, -end => $stop, -strand => $strand, ));
}
$bases= $dnaseq->subseq($sloc); # should always return forward strand dna, but for mixed exons
if($topstrand < 0 && $bases) { $bases= _dna_revcomp($bases); }
# if($topstrand < 0 && $bases) { warn "DEBUG: need _dna_revcomp(bases) for $name, loc=$baseloc; sloc=",$sloc->to_FTstring(),"\n"; }
$sloc->strand($topstrand); # BUT add in strand for to_FTstring ..
$gotloc= $sloc->to_FTstring(); # moved down; might mess w/ exon strands
## warn "feat=$name, got featloc=",$gotloc,"\n" if $DEBUG;
}
else {
## given warnings do we need to swap start,stop for -strand ?
my $sloc= new Bio::Location::Simple(
-start => $start, -end => $stop, -strand => $strand, );
$bases= $dnaseq->subseq($sloc);
$gotloc= $sloc->to_FTstring();
}
}
print STDERR "dna-file: $name, bases=",length($bases),"\n" if $DEBUG > 1;
if (!$subrange && ($gotloc ne $baseloc)) {
my $ok= 0;
## check for silly 123..123 => 123 change
while ($baseloc =~ m/(\d+)\.\.(\d+)/g) {
my($a,$b)= ($1,$2);
if ($a eq $b) { $ok=1; last; }
}
warn "dna-file: WARNING $name, loc-out=$gotloc ne loc-in=$baseloc\n"
unless $ok;
}
return $bases;
}
sub _dna_revcomp {
my($bases)= @_;
## Bioperl: $bases = Bio::PrimarySeq->new(-seq => $bases)->revcom()->seq();
## simple way:
$bases =~ tr/acgtrymkswhbvdnxACGTRYMKSWHBVDNX/tgcayrkmswdvbhnxTGCAYRKMSWDVBHNX/; # complement, doesnt do RNA u/U>t/T
return CORE::reverse $bases;
}
sub _isold {
my($source,$target) = @_;
## not for symlinks or dirs
my $res= 0;
my $targtime= -M $target; ## -M is file age in days.hrs before now
if (! -f $target) { return 1; }
elsif ( -l $source ) {
# $source= _getLinkOriginal($source);
$res= (-M $source) < $targtime;
}
elsif ( -f $source ) {
$res= (-M $source) < $targtime;
}
else { $res= 0; }
return $res;
}
sub _getLinkOriginal {
my($source) = @_;
my $rsource= readlink($source);
return $source unless ($rsource);
if ($rsource =~ m/^\.\./) {
my $at= rindex( $source,'/');
$rsource= substr($source,0,$at) . '/' . $rsource;
}
return $rsource;
}
## option to symlink ReleaseDir to 'current' for genomeweb path
sub makeCurrentLink
{
my ($self)= @_;
my $sok= 0;
my $reldir = $self->getReleaseDir();
my $subdir = $self->config->{relfull} || $self->config->{rel} || "release";
my $curdir = $reldir;
unless($curdir =~ s,$subdir$,current,) {
$sok= -1; return $sok; # what?
}
if( -l $curdir) {
my $lsource = _getLinkOriginal($curdir);
if ($lsource =~ m/$subdir$/) { $sok= 1; }
else {
warn "unlink old $curdir -> $lsource\n" if $DEBUG;
unlink($curdir); $sok= 1; # unlink err?
}
}
elsif( -d $curdir) {
warn "'$curdir' is directory, not symlink\n";
$sok= -1;
}
unless(-e $curdir) {
(my $topdir= $curdir) =~ s,current$,,;
my $olddir= $ENV{'PWD'}; #?? not safe?
$sok = chdir($topdir); # NOTE: relative link; subdir has no path
$sok = symlink($subdir,'current') if $sok;
chdir($olddir);
}
warn "Changed 'current' release symlink to $reldir; ok=$sok\n";
return $sok;
}
sub makePath
{
my($self, $dir, $errmsg, $failonerr)= @_;
return 0 unless($dir); # ok or not?
return 1 if(-d $dir);
eval { mkpath($dir,$DEBUG); }; # 0777 permission
if ($@) {
$failonerr= ($self->{failonerror}||$self->{automake}) unless(defined $failonerr);
warn "ERROR: Couldn't create path $dir: $@\n$errmsg\n";
die if ($failonerr);
return 0;
}
else { return 1; }
}
sub getReleaseDir
{
my($self)= @_;
my $config = $self->config;
my $releasedir= $config->{releasedir};
return $releasedir if ($releasedir && -d $releasedir);
## optiona/default: add full species to path, if not there ...
my $datadir= $config->{datadir}; # must exist
my $species= $self->speciesFull();
my $org= $self->config->{org};
my $subdir = $config->{relfull} || $config->{rel} || "release";
$datadir= catdir($datadir, $species) unless($datadir =~ m/$species|$org/);
$releasedir= catdir($datadir, $subdir);
$config->{releasedir} = $releasedir;
if( ! -d $datadir) { warn " missing data dir $datadir\n"; }
else { $self->makePath($releasedir, "Release path $subdir needed"); }
return $releasedir;
}
sub getReleaseSubdir
{
my($self, $subdir, $flags)= @_;
my $config= $self->config;
$flags ||= "";
unless(-d $subdir) {
my ($filename,$ext);
if ($subdir !~ m,/$, && $subdir =~ m,/, && $subdir =~ m,\.,) {
($filename, $subdir, $ext) = File::Basename::fileparse($subdir, '\.[^\.]+');
}
my $reldir= $self->getReleaseDir();
$subdir= catdir($reldir,$subdir) unless(-d $subdir);
$self->makePath($subdir, "Release subdir $subdir needed")
unless(-d $subdir || $flags =~ /nocreate|nomake/); ## mkpath
}
return $subdir;
}
# ## promote all to top of config ..
#
sub promoteRelease
{
my($self, $config)= @_;
my $cdate= $self->{date}; $cdate =~ s/[^\w-]/_/g; ## make sure no space in date !
my $relid= $config->{relid};
unless($relid) { $relid= $config->{relid}= $cdate; }
unless( ref $config->{release}->{$relid} ) {
$config->{release}->{$relid}= {
date => $cdate,
};
}
if($relid && ref $config->{release}->{$relid})
{
my %relh= %{$config->{release}->{$relid}};
foreach my $k (keys %relh) {
$config->{'release_'.$k}= $config->{$k}= $relh{$k} ;
# double-store: too many 'date' keys in conf files; use '${release_key}' by preference
}
}
unless($config->{rel}) { # used much, relid ~= rel
my $org= $config->{org} || "rel";
$config->{rel}= $org.$relid;
}
$config->{'release_id'}= $config->{rel} unless($config->{'release_id'}); ## alias of config variables ??
$self->updateConfigVars(); # install ${release_} variables
}
sub initData
{
my($self, $config, $oroot)= @_;
# check $self for params
unless(ref $config) { $config= $self->config || {}; }
$self->{config}= $config;
$self->{verbose}= $self->{verbose} || $config->{verbose};
## added $self->config->{newuser} ; turn on verbose, valid checks with this
$config->{no_csomesplit}= $self->{no_csomesplit} if($self->{no_csomesplit});
if (ref $config->{ENV}) {
foreach my $key (%{$config->{ENV}}) {
$ENV{$key}= $config->{ENV}->{$key} unless($ENV{$key});
}
}
$self->promoteRelease($config);
unless(defined $oroot && -d $oroot) {
if (defined $config->{ROOT}) { $oroot= $config->{ROOT}; }
elsif ($ENV{ARGOS_SERVICE_ROOT}) { $oroot= $ENV{ARGOS_SERVICE_ROOT}; }
elsif ($ENV{ARGOS_ROOT} && $config->{SERVICE}) { $oroot= $ENV{ARGOS_ROOT}.'/'.$config->{SERVICE}; }
elsif ($ENV{GMOD_ROOT}) { $oroot= $ENV{GMOD_ROOT}; }
unless(defined $oroot && -d $oroot) {
my $bin = "$FindBin::RealBin";
if ( -e "$bin/../common/") { $oroot= "$bin/../"; }
elsif ( -e "$bin/../conf/") { $oroot= "$bin/../"; }
# ^^ this is putting data into GMODTools/ folder - ok? no?
else { $oroot= "./"; }
$oroot=`cd "$oroot" && pwd`; chomp($oroot);
}
}
print STDERR "Using rootpath=$oroot\n" if $DEBUG; # is this bad?
$self->{rootpath} = $config->{rootpath} = $oroot; # gmod_root ??
my $datadir= $config->{datadir} || "genomes";
$datadir= "$oroot/$datadir" unless(-d $datadir);
if (!-d $datadir && -d $oroot) {
$self->makePath($datadir,
"** Need writeable data dir=$datadir\nChange configuration datadir\n", 1); #mkpath($datadir,$DEBUG);
}
$config->{datadir} = $datadir;
my $tmpdir= $config->{TMP} || $self->getReleaseSubdir( "tmp/"); # will make dirs inside datadir
unless( $tmpdir && $self->makePath($tmpdir,"** Need writeable TMP folder\nTrying other..\n",0) )
{
$tmpdir = File::Temp::tempdir( "gmodXXXX", TMPDIR => 1, CLEANUP => 1 );
$self->makePath($tmpdir,
"** Need writeable TMP=$tmpdir\nChange configuration TMP\n", 1);
warn "Using TMP=$tmpdir\n";
}
$config->{TMP} = $tmpdir;
$self->{idpattern}= $config->{idpattern} || '[A-Za-z]+\d+';
$fndel= $config->{filepart_delimiter} || '-';
my $fset= $self->getFilesetInfo('dna');
my $dnadir= $self->getReleaseSubdir( $fset->{path} || 'dna/');
$self->{dnadir} = $dnadir;
# see getChromosomeTable: $chromosome= $config->{chromosome} if (ref $config->{chromosome});
$self->{addids}= $config->{addids} if defined $config->{addids} ;
$self->{ignoredbresidues}= $config->{ignoredbresidues}
if defined $config->{ignoredbresidues} ;
## FIXME -- tests w/ this allfeats can be bad ...
@allfeats= (ref $config->{allfeats}) ? @{$config->{allfeats}}
: qw(
BAC CDS DNA_motif EST RNA_motif aberration_junction cDNA_clone enhancer five_prime_UTR
gene insertion_site intron mRNA mRNA_genscan mRNA_piecegenie mature_peptide ncRNA
oligo oligonucleotide point_mutation polyA_site processed_transcript protein protein_binding_site
pseudogene rRNA region regulatory_region repeat_region rescue_fragment
segment golden_path_fragment scaffold golden_path sequence_variant signal_peptide snRNA snoRNA
so source tRNA tRNA_trnascan three_prime_UTR transcription_start_site
transposable_element transposable_element_insertion transposable_element_insertion_site transposable_element_pred
);
%mapchr_pattern= %{ $config->{'mapchr_pattern'} } if ref $config->{'mapchr_pattern'};
# add all featset?
if (ref $config->{featset}) { @featset= @{$config->{featset}}; }
elsif ($config->{featset}) { @featset= ($config->{featset}); } # singleton
else { @featset= qw(gene mRNA CDS transcript translation
tRNA miscRNA transposon pseudogene gene_extended2000
five_prime_UTR three_prime_UTR intron
);
}
$config->{featset}= \@featset;
my @fastafeatok=();
push(@fastafeatok, @featset); # ?? not @allfeats
if (ref $config->{featmap}) {
my $fm= $config->{featmap};
foreach my $fk (keys %$fm) {
push(@fastafeatok, $fk);
if (ref $fm->{$fk} && defined $fm->{$fk}->{types}) {
my @ft= split(/[\s,;]/, $fm->{$fk}->{types} );
push(@fastafeatok, @ft);
}
}
}
$config->{fastafeatok}= \@fastafeatok;
## add these to %ENV before more configs so they get replaced ..
foreach my $k (@ENV_KEYS) { defined $$config{$k} and $ENV{$k}= $$config{$k}; }
}
#-----------
=head1
package Bio::GMOD::Bulkfiles::MySplitLocation
patch for Bio::Location::Split
-- moved to sep. file
=head1
package Bio::GMOD::Bulkfiles::MyLargePrimarySeq
-- moved to sep. file
patch to use Bio::Seq::LargePrimarySeq to read
feature locations from dna.raw files.
my $dnaseq= Bio::GMOD::Bulkfiles::MyLargePrimarySeq->new( -file => $dnafile);
$loc= new Bio::Location::something(...);
$bases= $dnaseq->subseq($loc);
=cut
1;