Raw content of Bio::SeqIO::locuslink
<![CDATA[# $Id: locuslink.pm,v 1.2.2.2 2003/03/13 02:09:20 lapp Exp $
#
# BioPerl module for Bio::SeqIO::locuslink
#
# Cared for by Keith Ching <kching at gnf.org>
#
# Copyright Keith Ching
#
# You may distribute this module under the same terms as perl itself
#
# (c) Keith Ching, kching at gnf.org, 2002.
# (c) GNF, Genomics Institute of the Novartis Research Foundation, 2002.
#
# You may distribute this module under the same terms as perl itself.
# Refer to the Perl Artistic License (see the license accompanying this
# software package, or see http://www.perl.com/language/misc/Artistic.html)
# for the terms under which you may use, modify, and redistribute this module.
#
# THIS PACKAGE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS OR IMPLIED
# WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
# MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
#
# POD documentation - main docs before the code
=head1 NAME
Bio::SeqIO::locuslink - DESCRIPTION of Object
=head1 SYNOPSIS
# don't instantiate directly - instead do
my $seqio = Bio::SeqIO->new(-format => "locuslink", -file => \STDIN);
=head1 DESCRIPTION
This module parses LocusLink into Bio::SeqI objects with rich
annotation, but no sequence.
The input file has to be in the LL_tmpl format - the tabular format
will not work.
The way the current implementation populates the object is rather a
draft work than a finished work of art. Note that at this stage the
locuslink entries cannot be round-tripped, because the parser loses
certain information. For instance, most of the alternative transcript
descriptions are not retained. The parser also misses any element
that deals with visual representation (e.g., 'button') except for the
URLs. Almost all of the pieces of the annotation are kept in the
L<Bio::Annotation::Collection> object.
=head1 FEEDBACK
=head2 Mailing Lists
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list. Your participation is much appreciated.
bioperl-l@bioperl.org - General discussion
http://bioperl.org/MailList.shtml - About the mailing lists
=head2 Reporting Bugs
Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via
email or the web:
bioperl-bugs@bioperl.org
http://bugzilla.bioperl.org/
=head1 AUTHOR - Keith Ching
Email kching at gnf.org
Describe contact details here
=head1 CONTRIBUTORS
Hilmar Lapp, hlapp at gmx.net
=head1 APPENDIX
The rest of the documentation details each of the object methods.
Internal methods are usually preceded with a _
=cut
package Bio::SeqIO::locuslink;
use strict;
use vars qw(@ISA);
use Bio::SeqIO;
use Bio::Seq::SeqFactory;
use Bio::Species;
use Bio::Annotation::DBLink;
#use Bio::Annotation::Reference;
use Bio::Annotation::Comment;
use Bio::Annotation::SimpleValue;
use Bio::Annotation::OntologyTerm;
use Bio::Annotation::Collection;
@ISA = qw(Bio::SeqIO);
# list of all the field names in locuslink
my @locuslink_keys = qw(
ACCNUM
ALIAS_PROT
ALIAS_SYMBOL
ASSEMBLY
BUTTON
CDD
CHR
COMP
CONTIG
CURRENT_LOCUSID
DB_DESCR
DB_LINK
ECNUM
EVID
EXTANNOT
GO
GRIF
LINK
LOCUSID
LOCUS_CONFIRMED
LOCUS_TYPE
MAP
MAPLINK
NC
NG
NM
NP
NR
OFFICIAL_GENE_NAME
OFFICIAL_SYMBOL
OMIM
ORGANISM
PHENOTYPE
PHENOTYPE_ID
PMID
PREFERRED_GENE_NAME
PREFERRED_PRODUCT
PREFERRED_SYMBOL
PRODUCT
PROT
RELL
STATUS
STS
SUMFUNC
SUMMARY
TRANSVAR
TYPE
UNIGENE
XG
XM
XP
XR
);
# list of fields to make simple annotations from
# fields not listed here or as a key in feature hash are ignored (lost).
my %anntype_map = (
SimpleValue => [qw(
ALIAS_PROT
ALIAS_SYMBOL
CDD
CHR
CURRENT_LOCUSID
ECNUM
EXTANNOT
MAP
NC
NR
OFFICIAL_GENE_NAME
OFFICIAL_SYMBOL
PHENOTYPE
PREFERRED_GENE_NAME
PREFERRED_PRODUCT
PREFERRED_SYMBOL
PRODUCT
RELL
SUMFUNC
)
],
Comment => [qw(
SUMMARY
)
],
);
# certain fields are not named the same as the symgene database list
my %dbname_map = (
pfam => 'Pfam',
smart => 'SMART',
NM => 'RefSeq',
NP => 'RefSeq',
XP => 'RefSeq',
XM => 'RefSeq',
NG => 'RefSeq',
XG => 'RefSeq',
XR => 'RefSeq',
PROT => 'GenBank',
ACCNUM => 'GenBank',
CONTIG => 'GenBank',
# certain fields are not named the same as the symgene
# database list: rename the fields the symgene database name
# key = field name in locuslink
# value = database name in sym
#GO => 'GO',
OMIM => 'MIM',
GRIF => 'GRIF',
STS => 'STS',
UNIGENE => 'UniGene',
);
# certain CDD entries use the wrong prefix for the accession number
# cddprefix will replace the key w/ the value for these entries
my %cddprefix = (
pfam => 'PF',
smart => 'SM',
);
# alternate mappings if one field does not exist
my %alternate_map = (
OFFICIAL_GENE_NAME => 'PREFERRED_GENE_NAME',
OFFICIAL_SYMBOL => 'PREFERRED_SYMBOL',
);
# for these field names, we only care about the first value X in value X|Y|Z
my @ll_firstelements = qw(
NM
NP
NG
XG
XM
XP
XR
PROT
STS
ACCNUM
CONTIG
GRIF
);
# these fields need to be flattened into a single string, using the given
# join string
my %flatten_tags = (
ASSEMBLY => '',
ORGANISM => '',
OFFICIAL_SYMBOL => '',
OFFICIAL_GENE_NAME => '',
LOCUSID => '',
PMID => '',
PREFERRED_SYMBOL => ', ',
PREFERRED_GENE_NAME => ', '
);
# set the default search pattern for all the field names
my %feature_pat_map = map { ($_ , "^$_: (.+)\n"); } @locuslink_keys;
sub _initialize {
my($self,@args) = @_;
$self->SUPER::_initialize(@args);
# overwrite the search pattern w/ the first value pattern
foreach my $key(@ll_firstelements){
$feature_pat_map{$key}="^$key: ([^|]+)";
}
# special search pattern for cdd entries
foreach my $key(keys %cddprefix) {
$feature_pat_map{$key}='^CDD: .+\|'.$key.'(\d+)';
}
# special patterns for specific fields
$feature_pat_map{MAP} = '^MAP: (.+?)\|';
$feature_pat_map{MAPHTML} = '^MAP: .+\|(<.+>)\|';
$feature_pat_map{GO} = '^GO: .+\|.+\|\w+\|(GO:\d+)\|';
$feature_pat_map{GO_DESC} = '^GO: .+\|(.+)\|\w+\|GO:\d+\|';
$feature_pat_map{GO_CAT} = '^GO: (.+)\|.+\|\w+\|GO:\d+\|';
$feature_pat_map{EXTANNOT} = '^EXTANNOT: (.+)\|(.+)\|\w+\|.+\|\d+';
# set the sequence factory of none has been set already
if(! $self->sequence_factory()) {
$self->sequence_factory(Bio::Seq::SeqFactory->new(
-type => 'Bio::Seq::RichSeq'));
}
}
#########################
#
sub search_pattern{
#
#########################
my ($self,
$entry, #text to search
$searchconfirm, #to make sure you got the right thing
$searchpattern,
$searchtype) = @_;
my @query = $entry=~/$searchpattern/gm;
if ($searchconfirm ne "FALSE"){
$self->warn("No $searchtype found\n$entry\n") unless @query;
foreach (@query){
if (!($_=~/$searchconfirm/)){
$self->throw("error\n$entry\n$searchtype parse $_ does not match $searchconfirm\n");
}
}#endforeach
}#endsearchconfirm
return(@query);
}#endsub
############
#
sub read_species{
#
############
my ($spline)=@_;
my $species;
my $genus;
($genus,$species)=$spline=~/([^ ]+) ([^ ]+)/;
my $make = Bio::Species->new();
$make->classification( ($species,$genus) );
return $make;
}
################
#
sub read_dblink{
#
################
my ($ann,$db,$ref)=@_;
my @results=$ref ? @$ref : ();
foreach my $id(@results){
if($id){
$ann->add_Annotation('dblink',
Bio::Annotation::DBLink->new(
-database =>$db ,
-primary_id =>$id));
}
}
return($ann);
}
################
#
sub read_reference{
#
################
my ($ann,$db,$results)=@_;
if($results){
chomp($results);
my @ids=split(/,/,$results);
$ann = read_dblink($ann,$db,\@ids) if @ids;
}
return $ann;
}#endsub
################
#
sub add_annotation{
#
################
my ($ac,$type,$text,$anntype)=@_;
my @args;
$anntype = 'SimpleValue' unless $anntype;
SWITCH : {
$anntype eq 'SimpleValue' && do {
push(@args, -value => $text, -tagname => $type);
last SWITCH;
};
$anntype eq 'Comment' && do {
push(@args, -text => $text, -tagname => 'comment');
last SWITCH;
};
}
$ac->add_Annotation("Bio::Annotation::$anntype"->new(@args));
return($ac);
}#endsub
################
#
sub add_annotation_ref{
#
################
my ($ann,$type,$textref)=@_;
my @text=$textref ? @$textref : ();
foreach my $text(@text){
$ann->add_Annotation($type,Bio::Annotation::SimpleValue->new(-value => $text));
}
return($ann);
}#endsub
################
#
sub make_unique{
#
##############
my ($ann,$key) = @_;
my %seen = ();
foreach my $dbl ($ann->remove_Annotations($key)) {
if(! $seen{$dbl->as_text()}) {
$seen{$dbl->as_text()} = 1;
$ann->add_Annotation($dbl);
}
}
return $ann;
}
################
#
sub next_seq{
#
##############
my ($self, @args)=@_;
my (@results,$search,$ref,$cddref);
# LOCUSLINK entries begin w/ >>
local $/="\n>>";
# slurp in a whole entry and return if no more entries
return unless my $entry = $self->_readline;
# strip the leading '>>' is it's the first entry
if (index($entry,'>>') == 0) { #first entry
$entry = substr($entry,2);
}
# we aren't interested in obsoleted entries, so we need to loop
# and skip those until we've found the next not obsoleted
my %record = ();
while($entry && ($entry =~ /\w/)) {
if (!($entry=~/LOCUSID/)){
$self->throw("No LOCUSID in first line of record. ".
"Not LocusLink in my book.");
}
# see whether it's an obsoleted entry, and if so jump to the next
# one entry right away
if($entry =~ /^CURRENT_LOCUSID:/m) {
# read next entry and continue
$entry = $self->_readline;
%record = ();
next;
}
# loop through list of features and get field values
# place into record hash as array refs
foreach my $key (keys %feature_pat_map){
$search=$feature_pat_map{$key};
@results=$self->search_pattern($entry,'FALSE',$search,$search);
$record{$key} = @results ? [@results] : undef;
}#endfor
# terminate loop as this one hasn't been obsoleted
last;
}
# we have reached the end-of-file ...
return unless %record;
# special processing for CDD entries like pfam and smart
my ($PRESENT,@keep);
foreach my $key(keys %cddprefix){
#print "check CDD $key\n";
if($record{$key}) {
@keep=();
foreach my $list (@{$record{$key}}) {
# replace AC with correct AC number
push(@keep,$cddprefix{$key}.$list);
}
# replace CDD ref with correctly prefixed AC number
$record{$key} = [@keep];
}
}
# modify CDD references @=();
if($record{CDD}) {
@keep=();
foreach my $cdd (@{$record{CDD}}) {
$PRESENT = undef;
foreach my $key (keys %cddprefix) {
if ($cdd=~/$key/){
$PRESENT = 1;
last;
}
}
push(@keep,$cdd) if(! $PRESENT);
}
$record{CDD} = [@keep];
}
# create annotation collection - we'll need it now
my $ann = Bio::Annotation::Collection->new();
foreach my $field(keys %dbname_map){
$ann=read_dblink($ann,$dbname_map{$field},$record{$field});
}
# add GO link as an OntologyTerm annotation
if($record{GO}) {
for(my $j = 0; $j < @{$record{GO}}; $j++) {
my $goann = Bio::Annotation::OntologyTerm->new(
-identifier => $record{GO}->[$j],
-name => $record{GO_DESC}->[$j],
-ontology => $record{GO_CAT}->[$j]);
$ann->add_Annotation($goann);
}
}
$ann=add_annotation_ref($ann,'URL',$record{LINK});
$ann=add_annotation_ref($ann,'URL',$record{DB_LINK});
# presently we can't store types of dblinks - hence make unique
make_unique($ann,'dblink');
# everything else gets a simple tag or comment value annotation
foreach my $anntype (keys %anntype_map) {
foreach my $key (@{$anntype_map{$anntype}}){
if($record{$key}){
foreach (@{$record{$key}}){
#print "$key\t\t$_\n";
$ann=add_annotation($ann,$key,$_,$anntype);
}
}
}
}
# flatten designated attributes into a scalar value
foreach my $field (keys %flatten_tags) {
if($record{$field}) {
$record{$field} = join($flatten_tags{$field},
@{$record{$field}});
}
}
# annotation that expects the array flattened out
$ann=read_reference($ann,'PUBMED',$record{PMID});
if($record{ASSEMBLY}) {
my @assembly=split(/,/,$record{ASSEMBLY});
$ann=read_dblink($ann,'GenBank',\@assembly);
}
# replace fields w/ alternate if original does not exist
foreach my $fieldval (keys %alternate_map){
if((! $record{$fieldval}) && ($record{$alternate_map{$fieldval}})){
$record{$fieldval}=$record{$alternate_map{$fieldval}};
}
}
# create sequence object (i.e., let seq.factory create one)
my $seq = $self->sequence_factory->create(
-verbose => $self->verbose(),
-accession_number => $record{LOCUSID},
-desc => $record{OFFICIAL_GENE_NAME},
-display_id => $record{OFFICIAL_SYMBOL},
-species => read_species($record{ORGANISM}),
-annotation => $ann);
# dump out object contents
# show_obj([$seq]);
return($seq);
}
################
#
sub show_obj{
#
################
my ($seqlistref)=@_;
my @list=@$seqlistref;
my $out = Bio::SeqIO->new('-fh' => \*STDOUT, -format => 'genbank' );
my ($ann,@values,$val);
foreach my $seq(@list){
$out->write_seq($seq);
$ann=$seq->annotation;
foreach my $key ( $ann->get_all_annotation_keys() ) {
@values = $ann->get_Annotations($key);
foreach my $value ( @values ) {
# value is an Bio::AnnotationI, and defines a "as_text" method
$val=$value->as_text;
print "Annotation ",$key,"\t\t",$val,"\n";
}
}
}
}#endsub
1;
]]>