my $self = shift;
my $slice = shift;
my $FH = shift;
my $FORMAT = shift;
#use only the core database to dump features (except for bloody snps)
my $lite = $slice->adaptor->db->remove_db_adaptor('lite');
my $meta = $slice->adaptor->db->get_MetaContainer;
#lump file handle and format string together for simpler method calls
my @ff = ($FH, $FORMAT);
my $value;
#source
my $classification = join(', ', @{$meta->get_classification()});
$self->write(@ff,'source', "1.." . $slice->length());
$self->write(@ff,'' , '/organism="'.$meta->get_scientific_name(). '"');
$self->write(@ff,'' , '/db_xref="taxon:'.$meta->get_taxonomy_id().'"');
#
# Transcripts & Genes
#
my @gene_slices;
if($self->is_enabled('gene')) {
push @gene_slices, $slice;
}
# Retrieve slices of other database where we need to pull genes from
my $gene_dbs = {'vegagene' => 'vega',
'estgene' => 'estgene'};
foreach my $gene_type (keys %$gene_dbs) {
if($self->is_enabled($gene_type)) {
my $db = $self->get_database($gene_dbs->{$gene_type});
if($db) {
my $sa = $db->get_SliceAdaptor();
push @gene_slices, $sa->fetch_by_name($slice->name());
} else {
warning("A [". $gene_dbs->{$gene_type} ."] database must be " .
"attached to this SeqDumper\n(via a call to " .
"attach_database) to retrieve genes of type [$gene_type]");
}
}
}
foreach my $gene_slice (@gene_slices) {
my @genes = @{$gene_slice->get_all_Genes(undef,undef, 1)};
while(my $gene = shift @genes) {
$value = $self->features2location( [$gene] );
$self->write( @ff, 'gene', $value );
$self->write( @ff, "", '/gene='.$gene->stable_id_version() );
if(defined($gene->display_xref)){
$self->write( @ff, "",'/locus_tag="'.$gene->display_xref->display_id.'"');
}
my $desc = $gene->description;
if(defined($desc) and $desc ne ""){
$desc =~ s/\"//;
$self->write( @ff, "", '/note="'.$gene->description.'"');
}
foreach my $transcript (@{$gene->get_all_Transcripts}) {
my $translation = $transcript->translation;
# normal transcripts get dumped differently than pseudogenes
if($translation) {
#normal transcript
$value = $self->features2location($transcript->get_all_Exons);
$self->write(@ff, 'mRNA', $value);
$self->write(@ff,'' , '/gene="'.$gene->stable_id_version().'"');
$self->write(@ff,''
,'/standard_name="'.$transcript->stable_id_version().'"');
# ...and a CDS section
$value =
$self->features2location($transcript->get_all_translateable_Exons);
$self->write(@ff,'CDS', $value);
my $codon_start = $self->transcript_to_codon_start($transcript);
$self->write(@ff,'' , qq{/codon_start="${codon_start}"}) if $codon_start > 1;
my $codon_table_id = $self->_get_codon_table_id($slice);
if($codon_table_id > 1){
$self->write(@ff,'' , '/transl_table='.$codon_table_id);
}
$self->write(@ff,'' , '/gene="'.$gene->stable_id_version().'"');
$self->write(@ff,'' , '/protein_id="'.$translation->stable_id_version().'"');
$self->write(@ff,'' ,'/note="transcript_id='.$transcript->stable_id_version().'"');
foreach my $dbl (@{$transcript->get_all_DBLinks}) {
my $db_xref = '/db_xref="'.$dbl->dbname().':'.$dbl->primary_id().'"';
my $go_db_xref = '/db_xref="'.$dbl->primary_id().'"';
$value = ($dbl->dbname()=~/GO/) ? $go_db_xref : $db_xref;
$self->write(@ff, '', $value);
}
my $pep = $transcript->translate();
if ($pep) {
$value = '/translation="'.$pep->seq().'"';
$self->write(@ff, '', $value);
}
} else {
#pseudogene
$value = $self->features2location($transcript->get_all_Exons);
$self->write(@ff, 'misc_RNA', $value);
$self->write(@ff,'' , '/gene="'.$gene->stable_id_version().'"');
foreach my $dbl (@{$transcript->get_all_DBLinks}) {
$value = '/db_xref="'.$dbl->dbname().':'.$dbl->primary_id().'"';
$self->write(@ff, '', $value);
}
$self->write(@ff,'' , '/note="'.$transcript->biotype().'"');
$self->write(@ff,''
,'/standard_name="'.$transcript->stable_id_version().'"');
}
}
}
# exons
foreach my $gene (@{$gene_slice->get_all_Genes(undef,undef,1)}) {
foreach my $exon (@{$gene->get_all_Exons}) {
$self->write(@ff,'exon', $self->features2location([$exon]));
$self->write(@ff,'' , '/note="exon_id='.$exon->stable_id_version().'"');
}
}
}
#
# genscans
#
if($self->is_enabled('genscan')) {
my @genscan_exons;
my @transcripts = @{$slice->get_all_PredictionTranscripts(undef,1)};
while(my $transcript = shift @transcripts) {
my $exons = $transcript->get_all_Exons();
push @genscan_exons, @$exons;
$self->write(@ff, 'mRNA', $self->features2location($exons));
$self->write(@ff, '', '/product="'.$transcript->translate()->seq().'"');
$self->write(@ff, '', '/note="identifier='.$transcript->stable_id_version.'"');
$self->write(@ff, '', '/note="Derived by automated computational' .
' analysis using gene prediction method:' .
$transcript->analysis->logic_name . '"');
}
}
#
# snps
#
if($self->is_enabled('variation') && $slice->can('get_all_VariationFeatures')) {
# $slice->adaptor->db->add_db_adaptor('lite', $lite) if $lite;
my @variations = @{$slice->get_all_VariationFeatures()};
while(my $snp = shift @variations) {
my $ss = $snp->start;
my $se = $snp->end;
#skip snps that hang off edge of slice
next if($ss < 1 || $se > $slice->length);
$self->write(@ff, 'variation', "$ss..$se");
$self->write(@ff, '' , '/replace="'.$snp->allele_string.'"');
#$self->write(@ff, '' , '/evidence="'.$snp->status.'"');
my $rs_id = $snp->variation_name();
my $db = $snp->source_name();
# foreach my $link ($snp->each_DBLink) {
# my $id = $link->primary_id;
# my $db = $link->database;
$self->write(@ff, '', "/db_xref=\"$db:$rs_id\"");
# }
}
# $slice->adaptor->db->remove_db_adaptor('lite') if $lite;
}
#
# similarity features
#
if($self->is_enabled('similarity')) {
foreach my $sim (@{$slice->get_all_SimilarityFeatures}) {
$self->write(@ff, 'misc_feature', $self->features2location([$sim]));
$self->write(@ff, '' , '/note="match: '.$sim->hseqname.
' : '.$sim->hstart.'..'.$sim->hend.'('.$sim->hstrand.')"');
}
}
#
# repeats
#
if($self->is_enabled('repeat')) {
my @rfs = @{$slice->get_all_RepeatFeatures()};
while(my $repeat = shift @rfs) {
$self->write(@ff, 'repeat_region', $self->features2location([$repeat]));
$self->write(@ff, '' , '/note="' . $repeat->repeat_consensus->name.
' repeat: matches ' . $repeat->hstart.'..'.$repeat->hend .
'('.$repeat->hstrand.') of consensus"');
}
}
#
# markers
#
if($self->is_enabled('marker') && $slice->can('get_all_MarkerFeatures')) {
my @markers = @{$slice->get_all_MarkerFeatures()};
while(my $mf = shift @markers) {
$self->write(@ff, 'STS', $self->features2location([$mf]));
if($mf->marker->display_MarkerSynonym) {
$self->write(@ff, '' , '/standard_name="' .
$mf->marker->display_MarkerSynonym->name . '"');
}
#grep out synonyms without a source
my @synonyms = @{$mf->marker->get_all_MarkerSynonyms};
@synonyms = grep {$_->source } @synonyms;
foreach my $synonym (@synonyms) {
$self->write(@ff, '', '/db_xref="'.$synonym->source.
':'.$synonym->name.'"');
}
$self->write(@ff, '', '/note="map_weight='.$mf->map_weight.'"');
}
}
#
# contigs
#
if($self->is_enabled('contig')) {
foreach my $segment (@{$slice->project('seqlevel')}) {
my ($start, $end, $slice) = @$segment;
$self->write(@ff, 'misc_feature',
$start .'..'. $end);
$self->write(@ff, '', '/note="contig '.$slice->seq_region_name .
' ' . $slice->start . '..' . $slice->end .
'(' . $slice->strand . ')"');
}
}
$slice->adaptor->db->add_db_adaptor('lite', $lite) if $lite;
}