Skip site navigation (1)Skip section navigation (2)

FreeBSD Manual Pages


home | help
Bio::SeqFeature::ToolsUserfContributed PBio::SeqFeature::Tools::Unflattener(3)

       Bio::SeqFeature::Tools::Unflattener - turns flat	list of
       genbank-sourced features	into a nested SeqFeatureI hierarchy

	 # standard / generic use - unflatten a	genbank	record
	 use Bio::SeqIO;
	 use Bio::SeqFeature::Tools::Unflattener;

	 # generate an Unflattener object
	 $unflattener =	Bio::SeqFeature::Tools::Unflattener->new;

	 # first fetch a genbank SeqI object
	 $seqio	=
	 my $out =
	 while ($seq = $seqio->next_seq()) {

	   # get top level unflattended	SeqFeatureI objects

	   @top_sfs = $seq->get_SeqFeatures;
	   foreach my $sf (@top_sfs) {
	       # do something with top-level features (eg genes)

       Most GenBank entries for	annotated genomic DNA contain a	flat list of
       features. These features	can be parsed into an equivalent flat list of
       Bio::SeqFeatureI	objects	using the standard Bio::SeqIO classes.
       However,	it is often desirable to unflatten this	list into something
       resembling actual gene models, in which genes, mRNAs and	CDSs are
       nested according	to the nature of the gene model.

       The BioPerl object model	allows us to store these kind of associations
       between SeqFeatures in containment hierarchies -- any SeqFeatureI
       object can contain nested SeqFeatureI objects. The
       Bio::SeqFeature::Tools::Unflattener object facilitates construction of
       these hierarchies from the underlying GenBank flat-feature-list

       For example, if you were	to look	at a typical GenBank DNA entry,	say,
       AE003644, you would see a flat list of features:


	 gene CG4491
	 mRNA CG4491-RA
	 CDS CG4491-PA

	 gene tRNA-Pro
	 tRNA tRNA-Pro

	 gene CG32954
	 mRNA CG32954-RA
	 mRNA CG32954-RC
	 mRNA CG32954-RB
	 CDS CG32954-PA
	 CDS CG32954-PB
	 CDS CG32954-PC

       These features have sequence locations, but it is not immediately clear
       how to write code such that each	mRNA is	linked to the appropriate CDS
       (other than relying on IDs which	is very	bad)

       We would	like to	convert	the above list into the	containment hierarchy,
       shown below:

	   mRNA	CG4491-RA
	     CDS CG4491-PA
	   tRNA	tRNA-Pro
	   mRNA	CG32954-RA
	     CDS CG32954-PA
	   mRNA	CG32954-RC
	     CDS CG32954-PC
	   mRNA	CG32954-RB
	     CDS CG32954-PB

       Where each feature is nested underneath its container. Note that	exons
       have been automatically inferred	(even for tRNA genes).

       We do this using	a call on a Bio::SeqFeature::Tools::Unflattener	object

	 @sfs =	$unflattener->unflatten_seq(-seq=>$seq);

       This would return a list	of the top level (i.e. container) SeqFeatureI
       objects - in this case, genes. Other top	level features are possible;
       for instance, the source	feature	which is always	present, and other
       features	such as	variation or misc_feature types.

       The containment hierarchy can be	accessed using the get_SeqFeature()
       call on any feature object - see	Bio::SeqFeature::FeatureHolderI.  The
       following code will traverse the	containment hierarchy for a feature:

	 sub traverse {
	   $sf = shift;	  #  $sf isa Bio::SeqfeatureI

	   # something with $sf!

	   # depth first traversal of containment tree
	   @contained_sfs = $sf->get_SeqFeatures;
	   traverse($_)	foreach	@contained_sfs;

       Once you	have built the hierarchy, you can do neat stuff	like turn the
       features	into 'rich' feature objects (eg
       Bio::SeqFeature::Gene::GeneStructure) or	convert	to a suitable format
       such as GFF3 or chadoxml	(after mapping to the Sequence Ontology); this
       step is not described here.

       Due to the quixotic nature of how features are stored in
       GenBank/EMBL/DDBJ, there	is no guarantee	that the default behaviour of
       this module will	produce	perfect	results. Sometimes it is hard or
       impossible to build a correct containment hierarchy if the information
       provided	is simply too lossy, as	is often the case. If you care deeply
       about your data,	you should always manually inspect the resulting
       containment hierarchy; you may have to customise	the algorithm for
       building	the hierarchy, or even manually	tweak the resulting hierarchy.
       This is explained in more detail	further	on in the document.

       However,	if you are satisfied with the default behaviour, then you do
       not need	to read	any further. Just make sure you	set the	parameter
       use_magic - this	will invoke incantations which will magically produce
       good results no matter what the idiosyncracies of the particular
       GenBank record in question.

       For example


       The success of this depends on the phase	of the moon at the time	the
       entry was submitted to GenBank. Note that the magical recipe is being
       constantly improved, so the results of invoking magic may vary
       depending on the	bioperl	release.

       If you are skeptical of magic, or you wish to exact fine	grained
       control over how	the entry is unflattened, or you simply	wish to
       understand more about how this crazy stuff works, then read on!

       Occasionally the	Unflattener will have problems with certain records.
       For example, the	record may contain inconsistent	data - maybe there is
       an exon entry that has no corresponding mRNA location.

       The default behaviour is	to throw an exception reporting	the problem,
       if the problem is relatively serious - for example, inconsistent	data.

       You can exert more fine grained control over this - perhaps you want
       the Unflattener to do the best it can, and report any problems. This
       can be done - refer to the methods.





       This is the default algorithm; you should be able to override any part
       of it to	customise.

       The core	of the algorithm is in two parts

       Partitioning the	flat feature list into groups
       Resolving the feature containment hierarchy for each group

       There are other optional	steps after the	completion of these two	steps,
       such as inferring exons;	we now describe	in more	detail what is going

   Partitioning	into groups
       First of	all the	flat feature list is partitioned into groups.

       The default way of doing	this is	to use the gene	attribute; if we look
       at two features from GenBank accession AE003644.3:

	    gene	    20111..23268
			    /note="last	curated	on Thu Dec 13 16:51:32 PST 2001"
	    mRNA	    join(20111..20584,20887..23268)

       Both these features share the same /gene	tag which is "noc", so they
       correspond to the same gene model (the CDS feature is not shown,	but
       this also has a tag-value /gene="noc").

       Not all groups need to correspond to gene models, but this is the most
       common use case;	later on we shall describe how to customise the

       Sometimes other tags have to be used; for instance, if you look at the
       entire record for AE003644.3 you	will see you actually need the use the
       /locus_tag attribute. This attribute is actually	not present in most

       You can override	this:

	 $collection->unflatten_seq(-seq=>$seq,	-group_tag=>'locus_tag');

       Alternatively, if you -use_magic, the object will try and make a	guess
       as to what the correct group_tag	should be.

       At the end of this step,	we should have a list of groups	- there	is no
       structure within	a group; the group just	serves to partition the	flat
       features. For the example data above, we	would have the following

	 [ source ]
	 [ gene	mRNA CDS ]
	 [ gene	mRNA CDS ]
	 [ gene	mRNA CDS ]
	 [ gene	mRNA mRNA mRNA CDS CDS CDS ]

       Multicopy Genes

       Multicopy genes are usually rRNAs or tRNAs that are duplicated across
       the genome. Because they	are functionally equivalent, and usually have
       the same	sequence, they usually have the	same group_tag (ie gene
       symbol);	they often have	a /note	tag giving copy	number.	This means
       they will end up	in the same group. This	is undesirable,	because	they
       are spatially disconnected.

       There is	another	step, which involves splitting spatially disconnected
       groups into distinct groups

       this would turn this

	[gene-rrn3 rRNA-rrn3 gene-rrn3 rRNA-rrn3]

       into this

	[gene-rrn3 rRNA-rrn3] [gene-rrn3 rRNA-rrn3]

       based on	the coordinates

       What next?

       The next	step is	to add some structure to each group, by	making
       containment hierarchies,	trees that represent how the features

   Resolving the containment hierarchy
       After the grouping is done, we end up with a list of groups which
       probably	contain	features of type 'gene', 'mRNA', 'CDS' and so on.

       Singleton groups	(eg the	'source' feature) are ignored at this stage.

       Each group is itself flat; we need to add an extra level	of
       organisation. Usually this is because different spliceforms
       (represented by the 'mRNA' feature) can give rise to different protein
       products	(indicated by the 'CDS'	feature). We want to correctly
       associate mRNAs to CDSs.

       We want to go from a group like this:

	 [ gene	mRNA mRNA mRNA CDS CDS CDS ]

       to a containment	hierarchy like this:


       In which	each CDS is nested underneath the correct corresponding	mRNA.

       For entries that	contain	no alternate splicing, this is simple; we know
       that the	group

	 [ gene	mRNA CDS ]

       Must resolve to the tree


       How can we do this in entries with alternate splicing? The bad news is
       that there is no	guaranteed way of doing	this correctly for any GenBank
       entry. Occasionally the submission will have been done in such a	way as
       to reconstruct the containment hierarchy. However, this is not
       consistent across databank entries, so no generic solution can be
       provided	by this	object.	This module does provide the framework within
       which you can customise a solution for the particular dataset you are
       interested in - see later.

       The good	news is	that there is an inference we can do that should
       produce pretty good results the vast majority of	the time. It uses
       splice coordinate data -	this is	the default behaviour of this module,
       and is described	in detail below.

   Using splice	site coordinates to infer containment
       If an mRNA is to	be the container for a CDS, then the splice site
       coordinates (or intron coordinates, depending on	how you	look at	it) of
       the CDS must fit	inside the splice site coordinates of the mRNA.

       Ambiguities can still arise, but	the results produced should still be
       reasonable and consistent at the	sequence level.	Look at	this fake

	 mRNA	 XXX---XX--XXXXXX--XXXX		join(1..3,7..8,11..16,19..23)
	 mRNA	 XXX-------XXXXXX--XXXX		join(1..3,11..16,19..23)
	 CDS		     XXXX--XX		join(13..16,19..20)
	 CDS		     XXXX--XX		join(13..16,19..20)

       [obviously the positions	have been scaled down]

       We cannot unambiguously match mRNA with CDS based on splice sites,
       since both CDS share the	splice site locations 16^17 and	18^19.
       However,	the consequences of making a wrong match are probably not very
       severe. Any annotation data attached to the first CDS is	probably
       identical to the	seconds	CDS, other than	identifiers.

       The default behaviour of	this module is to make an arbitrary call where
       it is ambiguous (the mapping will always	be bijective; i.e. one mRNA ->
       one CDS).

       [TODO: NOTE: not	tested on EMBL data, which may not be bijective; ie
       two mRNAs can share the same CDS??]

       This completes the building of the containment hierarchy; other
       optional	step follow

   Inferring exons from	mRNAs
       This step always	occurs if -use_magic is	invoked.

       In a typical GenBank entry, the exons are implicit. That	is they	can be
       inferred	from the mRNA location.

       For example:

	    mRNA	    join(20111..20584,20887..23268)

       This tells us that this particular transcript has two exons. In
       bioperl,	the mRNA feature will have a 'split location'.

       If we call


       This will generate the necessary	exon features, and nest	them under the
       appropriate mRNAs. Note that the	mRNAs will no longer have split
       locations - they	will have simple locations spanning the	extent of the
       exons. This is intentional, to avoid redundancy.

       Occasionally a GenBank entry will have both implicit exons (from	the
       mRNA location) and explicit exon	features.

       In this case, exons will	still be transferred. Tag-value	data from the
       explicit	exon will be transferred to the	implicit exon. If exons	are
       shared between mRNAs these will be represented by different objects.
       Any inconsistencies between implicit and	explicit will be reported.

       tRNAs and other noncoding RNAs

       exons will also be generated from these features

   Inferring mRNAs from	CDS
       Some GenBank entries represent gene models using	features of type gene,
       mRNA and	CDS; some entries just use gene	and CDS.

       If we only have gene and	CDS, then the containment hierarchies will
       look like this:


       If we want the containment hierarchies to be uniform, like this


       Then we must create an mRNA feature. This will have identical
       coordinates to the CDS. The assumption is that there is either no
       untranslated region, or it is unknown.

       To do this, we can call


       This is taken care of automatically, if -use_magic is invoked.

   Customising the grouping of features
       The default behaviour is	suited mostly to building models of protein
       coding genes and	noncoding genes	from genbank genomic DNA submissions.

       You can change the tag used to partition	the feature by passing in a
       different group_tag argument - see the unflatten_seq() method

       Other behaviour may be desirable. For example, even though SNPs
       (features of type 'variation' in	GenBank) are not actually part of the
       gene model, it may be desirable to group	SNPs that overlap or are
       nearby gene models.

       It should certainly be possible to extend this module to	do this.
       However,	I have yet to code this	part!!!	If anyone would	find this
       useful let me know.

       In the meantime,	you could write	your own grouping subroutine, and feed
       the results into	unflatten_groups() [see	the method documentation

   Customising the resolution of the containment hierarchy
       Once the	flat list of features has been partitioned into	groups,	the
       method unflatten_group()	is called on each group	to build a tree.

       The algorithm for doing this is described above;	ambiguities are
       resolved	by using splice	coordinates. As	discussed, this	can be

       Some submissions	may contain information	in tags/attributes that	hint
       as to the mapping that needs to be made between the features.

       For example, with the Drosophila	Melanogaster release 3 submission, we
       see that	CDS features in	alternately spliced mRNAs have a form like

	    CDS		    join(145588..145686,145752..146156,146227..146493)
			    /note="CG32954 gene	product	from transcript	CG32954-RA"

       Here the	/note tag provides the clue we need to link CDS	to mRNA
       (highlighted with ^^^^).	We just	need to	find the mRNA with the tag


       I have no idea how consistent this practice is across submissions; it
       is consistent for the fruitfly genome submission.

       We can customise	the behaviour of unflatten_group() by providing	our
       own resolver method. This obviously requires a bit of extra
       programming, but	there is no way	to get around this.

       Here is an example of how to pass in your own resolver; this example
       basically checks	the parent (container) /product	tag to see if it
       matches the required string in the child	(contained) /note tag.

					-resolver_method=>sub {
					    my $self = shift;
					    my ($sf, @candidate_container_sfs) = @_;
					    if ($sf->has_tag('note')) {
						my @notes = $sf->get_tag_values('note');
						my @trnames = map {/from transcript\s+(.*)/;
								   $1} @notes;
						@trnames = grep	{$_} @trnames;
						my $trname;
						if (@trnames ==	0) {
						elsif (@trnames	== 1) {
						    $trname = $trnames[0];
						else {
						    $self->throw("AMBIGUOUS: @trnames");
						my @container_sfs =
						  grep {
						      my ($product) =
							$_->has_tag('product') ?
							  $_->get_tag_values('product')	:
						      $product eq $trname;
						  } @candidate_container_sfs;
						if (@container_sfs == 0) {
						elsif (@container_sfs == 1) {
						    # we got it!
						    return $container_sfs[0];
						else {

       the resolver method is only called when there is	more than one

   Parsing mRNA	records
       Some of the entries in sequence databanks are for mRNA sequences	as
       well as genomic DNA. We may want	to build models	from these too.


       Open question - what would these	look like?

       Ideally we would	like a way of combining	a mRNA record with the
       corresponding SeFeature entry from the appropriate genomic DNA record.
       This could be problemmatic in some cases	- for example, the mRNA
       sequences may not match 100% (due to differences	in strain, assembly
       problems, sequencing problems, etc). What then...?

       Feature table description

   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 lists  Your participation is much appreciated.			       - General discussion	- About	the mailing lists

       Please direct usage questions or	support	issues to the mailing list:

       rather than to the module maintainer directly. Many experienced and
       reponsive experts will be able look at the problem and quickly address
       it. Please include a thorough description of the	problem	with code and
       data examples if	at all possible.

   Reporting Bugs
       report bugs to the Bioperl bug tracking system to help us keep track
       the bugs	and their resolution.  Bug reports can be submitted via	the

AUTHOR - Chris Mungall

       The rest	of the documentation details each of the object	methods.
       Internal	methods	are usually preceded with a _

	Title	: new
	Usage	: $unflattener = Bio::SeqFeature::Tools::Unflattener->new();
	Function: constructor
	Example	:
	Returns	: a new	Bio::SeqFeature::Tools::Unflattener
	Args	: see below


	 -seq	    : A	L<Bio::SeqI> object (optional)
		      the sequence to unflatten; this can also be passed in
		      when we call unflatten_seq()

	 -group_tag : a	string representing the	/tag used to partition flat features
		      (see discussion above)

	Title	: seq
	Usage	: $unflattener->seq($newval)
	Example	:
	Returns	: value	of seq (a Bio::SeqI)
	Args	: on set, new value (a Bio::SeqI, optional)

       The Bio::SeqI object should hold	a flat list of Bio::SeqFeatureI
       objects;	this is	the list that will be unflattened.

       The sequence object can also be set when	we call	unflatten_seq()

	Title	: group_tag
	Usage	: $unflattener->group_tag($newval)
	Example	:
	Returns	: value	of group_tag (a	scalar)
	Args	: on set, new value (a scalar or undef,	optional)

       This is the tag that will be used to collect elements from the flat
       feature list into groups; for instance, if we look at two typical
       GenBank features:

	    gene	    20111..23268
			    /note="last	curated	on Thu Dec 13 16:51:32 PST 2001"
	    mRNA	    join(20111..20584,20887..23268)

       We can see that these comprise the same gene model because they share
       the same	/gene attribute; we want to collect these together in groups.

       Setting group_tag is optional. The default is to	use 'gene'. In the
       example above, we could also use	/locus_tag

	Title	: partonomy
	Usage	: $unflattener->partonomy({mRNA=>'gene', CDS=>'mRNA')
	Example	:
	Returns	: value	of partonomy (a	scalar)
	Args	: on set, new value (a scalar or undef,	optional)

       A hash representing the containment structure that the seq_feature
       nesting should conform to; each key represents the contained (child)
       type; each value	represents the container (parent) type.

	Title	: structure_type
	Usage	: $unflattener->structure_type($newval)
	Example	:
	Returns	: value	of structure_type (a scalar)
	Args	: on set, new value (an	int or undef, optional)

       GenBank entries conform to different flavours, or structure types. Some
       have mRNAs, some	do not.

       Right now there are only	two base structure types defined. If you set
       the structure type, then	appropriate unflattening action	will be	taken.
       The presence or absence of explicit exons does not affect the structure

       If you invoke -use_magic	then this will be set automatically, based on
       the content of the record.

       Type 0 (DEFAULT)
	   typically contains


	   with	this structure type, we	want the seq_features to be nested
	   like	this


	   exons and introns are implicit from the mRNA	'join' location

	   to get exons	from the mRNAs,	you will need this call	(see below)


       Type 1
	   typically contains

	     exon [optional]
	     intron [optional]

	   there are no	mRNA features

	   with	this structure type, we	want the seq_features to be nested
	   like	this


	   exon	and intron may or may not be present; they may be implicit
	   from	the CDS	'join' location

	Title	: get_problems
	Usage	: @probs = get_problems()
	Function: Get the list of problem(s) for this object.
	Example	:
	Returns	: An array of [severity, description] pairs
	Args	:

       In the course of	unflattening a record, problems	may occur. Some	of
       these problems are non-fatal, and can be	ignored.

       Problems	are represented	as arrayrefs containing	a pair [severity,

       severity	is a number, the higher, the more severe the problem

       the description is a text string

	Title	: clear_problems
	Usage	:
	Function: resets the problem list to empty
	Example	:
	Returns	:
	Args	:

	Title	: report_problems
	Usage	: $unflattener->report_problems(\*STDERR);
	Example	:
	Returns	:
	Args	: FileHandle (defaults to STDERR)

	Title	: ignore_problems
	Usage	: $obj->ignore_problems();
	Example	:
	Returns	:
	Args	:

       Unflattener is very particular about problems it	finds along the	way.
       If you have set the error_threshold such	that less severe problems do
       not cause exceptions, Unflattener still expects you to
       report_problems() at the	end, so	that the user of the module is aware
       of any inconsistencies or problems with the data. In fact, a warning
       will be produced	if there are unreported	problems. To silence, this
       warning,	call the ignore_problems() method before the Unflattener
       object is destroyed.

	Title	: error_threshold
	Usage	: $obj->error_threshold($severity)
	Example	:
	Returns	: value	of error_threshold (a scalar)
	Args	: on set, new value (an	integer)

       Sets the	threshold above	which errors cause this	module to throw	an
       exception. The default is 0; all	problems with a	severity > 0 will
       cause an	exception.

       If you raise the	threshold to 1,	then the unflattening process will be
       more lax; problems of severity==1 are generally non-fatal, but may
       indicate	that the results should	be inspected, for example, to make
       sure there is no	data loss.

	Title	: unflatten_seq
	Usage	: @sfs = $unflattener->unflatten_seq($seq);
	Function: turns	a flat list of features	into a list of holder features
	Example	:
	Returns	: list of Bio::SeqFeatureI objects
	Args	: see below

       partitions a list of features then arranges them	in a nested tree; see
       above for full explanation.

       note - the Bio::SeqI object passed in will be modified


	 -seq	:	   a Bio::SeqI object; must contain Bio::SeqFeatureI objects
			   (this is optional if	seq has	already	been set)

	 -use_magic:	   if TRUE (ie non-zero) then magic will be invoked;
			   see discussion above.

	 -resolver_method: a CODE reference
			   see the documentation above for an example of
			   a subroutine	that can be used to resolve hierarchies
			   within groups.

			   this	is optional - if nothing is supplied, a	default
			   subroutine will be used (see	below)

	 -group_tag:	   a string
			   [ see the group_tag() method	]
			   this	overrides the default group_tag	which is 'gene'

	Title	: unflatten_groups
	Usage	:
	Function: iterates over	groups,	calling	unflatten_group() [see below]
	Example	:
	Returns	: list of Bio::SeqFeatureI objects that	are holders
	Args	: see below


	 -groups:	   list	of list	references; inner list is of Bio::SeqFeatureI objects
			   e.g.	 ( [$sf1], [$sf2, $sf3,	$sf4], [$sf5, ...], ...)

	 -resolver_method: a CODE reference
			   see the documentation above for an example of
			   a subroutine	that can be used to resolve hierarchies
			   within groups.

			   this	is optional - a	default	subroutine will	be used

       NOTE: You should	not need to call this method, unless you want fine
       grained control over how	the unflattening process.

	Title	: unflatten_group
	Usage	:
	Function: nests	a group	of features into a feature containment hierarchy
	Example	:
	Returns	: Bio::SeqFeatureI objects that	holds other features
	Args	: see below


	 -group:	   reference to	list of	Bio::SeqFeatureI objects

	 -resolver_method: a CODE reference
			   see the documentation above for an example of
			   a subroutine	that can be used to resolve hierarchies
			   within groups

			   this	is optional - a	default	subroutine will	be used

       NOTE: You should	not need to call this method, unless you want fine
       grained control over how	the unflattening process.

	Title	: feature_from_splitloc
	Usage	: $unflattener->feature_from_splitloc(-features=>$sfs);
	Example	:
	Returns	:
	Args	: see below

       At this time all	this method does is generate exons for mRNA or other
       RNA features


	 -feature:    a	Bio::SeqFeatureI object	(that conforms to Bio::FeatureHolderI)
	 -seq:	      a	Bio::SeqI object that contains Bio::SeqFeatureI	objects
	 -features:   an arrayref of Bio::SeqFeatureI object

	Title	: infer_mRNA_from_CDS
	Usage	:
	Example	:
	Returns	:
	Args	:

       given a "type 1"	containment hierarchy


       this will infer the uniform "type 0" containment	hierarchy


       all the children	of the CDS will	be moved to the	mRNA

       a "type 2" containment hierarchy	is mixed type "0" and "1" (for
       example,	see

	Title	: remove_types
	Usage	: $unf->remove_types(-seq=>$seq, -types=>["mRNA"]);
	Example	:
	Returns	:
	Args	:

       removes features	of a set type

       useful for pre-filtering	a genbank record; eg to	get rid	of STSs

       also, there is no way to	unflatten UNLESS the bogus
       mRNAs in	these records are removed (or changed to a different type) -
       they just confuse things	too much

perl v5.32.0			  2019-1Bio::SeqFeature::Tools::Unflattener(3)


Want to link to this manual page? Use this URL:

home | help