[Bioperl-l] Re: GFF3 preliminary
Richard Durbin
rd at sanger.ac.uk
Thu Feb 20 16:26:41 EST 2003
Here are more general comments on the GFF 3 proposal:
1) I'd like to add the following comment to the definition of columns 4
and 5 (start and end):
When the feature is an interbase junction, such as a splice site or
insertion point, by convention the start and end should be the
flanking two bases.
2) Following the exchange triggered by Ian's email, I suggest the
definition for column 8 (phase) replaces "start of the feature" by
"5' end of the feature (start if forward strand, end if reverse
strand)". And the end of this section could point out that if no strand
is given phase is not meaningful and should be '.'.
3) I'm not convinced by the format for the Align string. This requires
a character per aligned base. There are a variety of run-length type
encodings in common use that are much more compact. e.g. Ensembl uses a
string such as "60M1D8M3I15M" to mean "60 match, then 1 delete, then 8
match, then 3 insert, then 15 match". They call this CIGAR, but when I
talked to Guy Slater, who invented CIGAR for exonerate, his version is
subtly different: "M 60 D 1 M 8 I 3 M 15" for the same string (see
http://www.ensembl.org/Docs/wiki/html/EnsemblDocs/CigarFormat.html).
Jim Kent also has something like this. I'd prefer us to standardise on
one of these formats, all of which are very short for ungapped matches.
Otherwise alignments will be too large. I don't think we need to know
about mismatches in the format. If there are regions unaligned in both
sequences we can simply have xIyD or the equivalent in whatever syntax
is agreed. We ought to agree how to handle split codons in DNA to
protein alignments as part of this spec.
4) After discussion with Michele last week I propose we don't introduce
relative coordinates with respect to an ID. First, this adds
significant complexity that you are requiring any parser to handle, even
one dealing with simple features not genes with exons and coding
sequences. I bet many programs that read gff would not handle it.
Second it creates the potential for getting things wrong, particularly
with IDs corresponding to objects on the reverse strand. Third, I don't
see that it saves effort - someone has to map backwards and forwards
from relative to absolute coordinates and I don't see that GFF, which we
all agree should be kept simple as a primary goal, should carry that
burden. It should be done elsewhere. Fourth, what is the scope? Can
you write -1..2000 for an ID that is 1000bp long? (If so this had better
mean 2 bases before the start of the object, but that might confuse some
people.) Fifth, it adds interdependency between lines, making actions
like subsetting using grep or simple perl much more likely to break.
5) I find the statement about disjunct coordinates being implied by
sharing the same ID at odds with the statement that the ID is unique
within the scope of the file. In particular what would happen if this
ID got used as the seqid for new features (although I argued above that
we should not go that way). I would much prefer that GFF features are
only one per line, and are all simple intervals or interbase junctions.
When we want to represent something more complex we should agree on
how to do it out of these primitives, not define GFF to itself represent
more complex structures.
Personally I would prefer to have one CDS line span the coding region,
with exons that (at least in part) contribute to it as children. e.g.
5 ctg123 flybase gene 43733 44677 . + .
ID=gene00001;Alias=ADAM1;Note=unc-3;GO_term=GO:12345,GO:33421
6 ctg123 flybase mRNA 43733 44677 . + .
ID=mRNA00001;Alias=ADAM1.t1;Parent=gene00001
7 ctg123 flybase mRNA 43733 44677 . + .
ID=mRNA00002;Alias=ADAM1.t2;Parent=gene00001
14 ctg123 flybase cds 43740 44677 . + 0
ID=cds00001;Parent=mRNA00001,gene00001
15 ctg123 flybase cds 43740 44677 . + 0
ID=cds00001;Parent=mRNA00002,gene00001
8 ctg123 flybase exon 43733 43961 . + .
ID=exon00001;Parent=mRNA00001,mRNA00002,cds00001,cds00002
9 ctg123 flybase exon 44030 44234 . + .
ID=exon00002;Parent=mRNA00001,mRNA00002,cds00001,cds00002
10 ctg123 flybase exon 44281 44328 . + .
ID=exon00003;Parent=mRNA00002,cds00002
11 ctg123 flybase exon 44521 44677 . + .
ID=exon00004;Parent=mRNA00001,mRNA00002,cds00001,cds00002
12 ctg123 flybase coding_start 43739 43740 . + .
Parent=cds00001,cds00002
13 ctg123 flybase coding_end 44677 44678 . + .
Parent=cds00001,cds00002
Note that I changed the coordinates of the coding_start and coding_end
junctions to flank the junction, and made them children of the CDS, not
the exon.
I prefer this because it allows one CDS to corresond to multiple mRNAs,
which is quite common with alternate promoters and poly-adenylation
sites. However I don't think it corresponds to Ensembl's view which
makes a CDS on exons within an mRNA, although it is compatible. I think
it does fit with Flybase's ideas, at least as once explained to me by
Bill Gelbart, and with what I would like to do for WormBase, though I
don't think we are all agreed on that yet. As for the genefinders, they
almost all just predict cds's and so it simplifies considerably because
there is no need for the mRNA features or identifiers - you just have
cds's as children of genes.
Note that this agrees with Lincoln's prescription at the end of his mail
that coding_start and _end information should not be with respect to the
exon.
Note also that by including coding_start and coding_end I was intending
to imply that there is an ATG and a stop codon encoded. i.e. if they
are missing then this might just be a gene fragment.
I know there are many ways to handle these things. We should try to
converge on a standard, but I think this may be an adjunct document, not
part of the GFFv3 standard itself, just as GTF was to GFFv2.
Richard
Lincoln Stein wrote:
> Hi,
>
> Following up on discussions with Jim Kent, Suzi Lewis, Michele Clamp
> and Richard Durbin, here is a new version of the GFF3 proposal.
>
> Suzi, could you post this to song.sourceforge.net, when you have a
> chance? I don't seem to have write permissions to the htdocs
> directory.
>
> Best,
>
> Lincoln
>
>
> GENERIC FEATURE FORMAT VERSION 3: A PROPOSAL
>
> Author: Lincoln Stein
> Date: 19 February 2003
> Version: 0.2
>
> Although there are many richer ways of representing genomic features
> via XML, the stubborn persistence of a variety of ad-hoc tab-delimited
> flat file formats declares the bioinformatics community's need for a
> simple format that can be modified with a text editor and processed
> with shell tools like grep. The GFF format, although widely used, has
> fragmented into multiple incompatible dialects. When asked why they
> have modified the published Sanger specification, bioinformaticists
> frequently answer that the format was insufficient for their needs,
> and they needed to extend it. The proposed GFF3 format addresses the
> most common extensions to GFF, while preserving backward compatibility
> with previous formats. The new format:
>
> 1) adds a mechanism for representing more than one level
> of hierarchical grouping of features and subfeatures.
> 2) separates the ideas of group membership and feature name/id
> 3) constrains the feature type field to be taken from a controlled
> vocabulary.
> 4) allows a single feature, such as an exon, to belong to more than
> one group at a time.
> 5) one level of relative addressing for subfeatures (e.g. exons
> can be expressed in transcript coordinates)
> 6) an explicit convention for pairwise alignments
> 7) an explicit convention for features that occupy disjunct regions
>
> The format consists of 10 columns, separated by spaces. The following
> unescaped characters are allowed within fields:
> [a-zA-Z0-9.:;=%^*$@!+_?-]. All other characters must must be escaped
> using the URL escaping conventions. Unescaped quotation marks,
> backslashes and other ad-hoc escaping conventions that have been added
> to the GFF format are explicitly forbidden. The =, ; and % characters
> have reserved meanings as described below.
>
> Undefined fields are replaced with the "." character, as described in
> the original GFF spec.
>
> Column 1: "seqid"
>
> The ID of the landmark used to establish the coordinate system for the
> current feature. IDs must contain alphanumeric characters.
> Whitespace, if present, must be escaped using URL escaping rules
> (e.g. space="%20" or "+").
>
> Column 2: "source"
>
> The source of the feature. This is unchanged from the older GFF specs
> and is not part of a controlled vocabulary.
>
> Column 3: "type"
>
> The type of the feature (previously called the "method"). This is
> constrained to be either: (a) a term from SOFA; or (b) a SOFA
> accession number. The latter alternative is distinguished using the
> syntax SOFA:000000.
>
> Columns 4 & 5: "start" and "end"
>
> The start and end of the feature, in 1-based integer coordinates,
> relative to the landmark given in column 1. Start is less than end.
>
> Column 6: "score"
>
> The score of the feature, a floating point number. As in earlier
> versions of the format, the semantics of the score are ill-defined.
> It is strongly recommended that E-values be used for sequence
> similarity features, and that P-values be used for ab initio gene
> prediction features.
>
> Column 7: "strand"
>
> The strand of the feature. + for positive strand (relative to the
> landmark), - for minus strand, and . for features that are not
> stranded. In addition, ? can be used for features whose strandedness
> is relevant, but unknown.
>
> Column 8: "phase"
>
> The phase of the feature, for protein-encoding featues (primarily
> CDSs). This is an integer-valued field with the values 0, 1, or 2.
> The integer indicates the offset from the start of the feature to the
> first base of the first codon in the reading frame. "." is used for
> features that do not corresponding to a reading frame.
>
> Column 9: "attributes"
>
> A list of feature attributes in the format tag=value. Multiple
> tag=value pairs are separated by semicolons. URL escaping rules are
> used for tags or values containing the following characters: ",=;".
> Whitespace should be replaced with the "+" character or the %20 URL
> escape. This will allow the file to survive text processing programs
> that convert tabs into spaces.
>
> Five tags are predefined:
>
> ID Indicates the name of the feature. IDs must be unique
> within the scope of the GFF file.
>
> Alias A descriptive name for the feature. It is suggested that
> this tag be used whenever a secondary identifier for the
> feature is needed, such as display names, locus names and
> accession numbers. Unlike ID, there is no requirement
> that Alias be unique within the file.
>
> Parent Indicates the parent of the feature. A parent ID can be
> used to group exons into transcripts, transcripts into
> genes, an so forth. A feature may have multiple parents.
>
> Target Indicates the target of a nucleotide to nucleotide or
> nucleotide to protein alignment. The format of the
> value is "target_id:start..end" Start may be greater
> than end to indicate a + strand alignment to the
> reverse complement of a target nucleotide sequence.
>
> Align The alignment of the feature to the target if the two
> are not colinear. The alignment is a string containing
> the four characters "|X^v", where "|" indicates an
> aligned match, "X" indicates an aligned mismatch, "^"
> indicates a gap in the feature, and "v" indicates a
> gap in the target.
>
> Multiple attributes of the same type are indicated by separating the
> values with the comma "," character, as in:
>
> Parent=AF2312,AB2812,abc-3
>
> Note that attribute names are case sensitive. "Parent" is not the
> same as "parent".
>
> In the example GFF3 file given below, the first column contains line
> numbers that I have added for the purposes of the narrative. Here are
> some common scenarios that I have attempted to illustrate:
>
> A) a simple feature, no public ID
>
> Line 2 in the example is a feature of type "repeat". It is located on
> the coordinate system defined by feature "ctg123", has a start and an
> end and no ID. It has an attribute named "Note" with value "ALU3."
>
> B) a simple feature with a public ID
>
> Line 3 is a feature of type clone. It has a start and an end. Its
> parent is undefined (no Parent attribute), but it has an ID attribute
> of "clone00001" and an Alias of "cTel33B."
>
> C) a feature with multiple attributes
>
> Line 5 is a feature of type "gene." It has no parent, and has
> attributes of type ID, Note, and GO_term.
>
> D) a hierarchical grouping of features
>
> Lines 5-13 demonstrate a hierarchical grouping. At the top level is
> line 5, which defines the extent of a "gene" with ID "gene00001".
> Below this are two features of type mRNA (lines 6 and 7). Their
> Parent attributes are set to "gene00001", indicating that this feature
> is their immediate parent. Their IDs are indicated as separate
> attributes.
>
> This pattern is repeated for the exons listed on lines 8-11. Exons
> exon00001, exon00002, and exon00004 belong to both of the transcripts.
> Therefore, their Parent attribute contains both the mRNA00001 and
> mRNA00002 IDs separated by a comma.
>
> Exon exon00003 belongs to mRNA00002 only, and therefore that
> transcript's ID is listed as the sole Parent.
>
> Lines 12 and 13 indicate coding_start and coding_end features. These
> subfeatures are hierarchically grouped underneath their corresponding
> exons, but they do not have independent public IDs.
>
> E) Disjunct coordinates
>
> Lines 14-16 illustrates a single feature -- the CDS corresponding to
> mRNA mRNA00001 -- which occupies multiple disjunct regions. The
> Parent attribute indicates that the CDS features belong to mRNA00001.
> However, the attribute column assigns each of lines 14-16 the same ID.
> Because the ID is the same, this is interpreted as a single feature
> that spans multiple disjunct coordinate ranges.
>
> NOTE: See "Representing Translations" for a discussion of why it might
> not be a good idea to use represent translations in this way.
>
> F) Alignments
>
> Lines 17-19 demonstrate an alignment of two sequences using the
> reserved Target attribute. Each non-gapped segment becomes a line in
> the GFF3 file. The segments each share the same ID, thereby
> indicating that the segments are disjunct regions of the same feature.
> The Target attribute indicates the ID of the target sequence (which
> does not have to be represented in the GFF3 file) and the start and
> end coordinates of the aligned target.
>
> Line 20 shows a gapped alignment using the Align attribute. This
> attribute's value should be interpreted this way:
>
>
> 1501 gatt*ctccc 1510 ctg123
> ||||^||X||
> 2001 gatttctgcc 2011 af923
>
> Unlike the GFF1 and GFF2 formats, the Parent attribute for gapped
> alignments may be empty. However, a valid alternative representation
> is to create a single "match" feature, and a series of "hsp" features
> contained within it. Lines 21-23 show this alternative
> representation.
>
> G) Relative coordinates
>
> Lines 24-27 illustrate using relative coordinate addressing in
> feature/subfeature relationships. Line 24 defines an mRNA that is
> positioned on sequence landmark "ctg123" from positions 5000 to 6000.
> Its ID field indicates that is mRNA03. Lines 25-27 are exon
> subfeatures of mRNA03 as indicated by their Parent attribute.
> However, the seqid field specifies mRNA03 as the parent coordinate
> system, thereby allowing the exons to begin at position 1.
>
> 0 ##gff-version 3
> 1 ##sequence-region ctg123:1..1497228
>
> 2 ctg123 flybase repeat 5000 5100 . . . Note=ALU3
> 3 ctg123 flybase clone 1 2679 . + . ID=clone00001;Alias=cTel33B
> 4 ctg123 flybase contig 1 1497228 . + . ID=contig0001;Alias=ctg123
>
> 5 ctg123 flybase gene 43733 44677 . + . ID=gene00001;Alias=ADAM1;Note=unc-3;GO_term=GO:12345,GO:33421
> 6 ctg123 flybase mRNA 43733 44677 . + . ID=mRNA00001;Alias=ADAM1.t1;Parent=gene00001
> 7 ctg123 flybase mRNA 43733 44677 . + . ID=mRNA00002;Alias=ADAM1.t2;Parent=gene00001
> 8 ctg123 flybase exon 43733 43961 . + . ID=exon00001;Parent=mRNA00001,mRNA00002
> 9 ctg123 flybase exon 44030 44234 . + . ID=exon00002;Parent=mRNA00001,mRNA00002
> 10 ctg123 flybase exon 44281 44328 . + . ID=exon00003;Parant=mRNA00002
> 11 ctg123 flybase exon 44521 44677 . + . ID=exon00004;Parent=mRNA00001,mRNA00002
> 12 ctg123 flybase coding_start 43740 43740 . + . Parent=exon00001
> 13 ctg123 flybase coding_end 44677 44677 . + . Parent=exon00004
>
> 14 ctg123 flybase cds 43740 43961 . + 0 ID=cds00001;Parent=mRNA00001
> 15 ctg123 flybase cds 44030 44234 . + 1 ID=cds00001;Parent=mRNA00001
> 16 ctg123 flybase cds 44521 44677 . + 1 ID=cds00001;Parent=mRNA00001
>
> 17 ctg123 flybase match 1 100 100 . . ID=match0001;Target=af923:1001..1100
> 18 ctg123 flybase match 101 500 80 . . ID=match0001;Target=af923:1101..1500
> 19 ctg123 flybase match 501 1000 80 . . ID=match0001;Target=af923:1501..2000
> 20 ctg123 flybase match 1501 1510 60 . . ID=match0001;Target=af923:2001..2011;Align=||||^||X||
>
> 21 ctg123 flybase match 5001 6000 100 . . ID=match0002;Target=ua388:1..1000
> 22 ctg123 flybase hsp 5001 5500 . . . Parent=match0002;Target=ua388:1..500
> 23 ctg123 flybase hsp 5501 6000 . . . Parent=match0002;Target-ua388:501.1000
>
> 24 ctg123 flybase mRNA 5000 6000 + . . ID=mRNA03;Alias=EVE1.t1
> 25 mRNA03 flybase exon 1 300 + . . ID=exon00005;Parent=mRNA03
> 26 mRNA03 flybase exon 301 400 + . . ID=exon00006;Parent=mRNA03
> 27 mRNA03 flybase exon 401 1000 + . . ID=exon00007;Parent=mRNA03
>
> =================================================================
>
> OTHER SYNTAX:
>
> Comments are preceded by the # symbol. Meta-data and directives are
> preceded by ##. The following directives are recognized:
>
> ##gff-version 3
> The GFF version, always 3 in this spec. This must
> be the topmost line of the file.
>
> ##sequence-region seqid:start..end
> The sequence segment referred to
> by this file, in the format seqid:start..end.
> This element is optional. If it occurs, it must be
> the second line of the file.
>
> ###
> This directive (three # signs in a row) indicates that all
> forward references to feature IDs that have been seen to this
> point have been resolved. After seeing this directive, a
> program that is processing the file serially can close off any
> open objects that it has created and return them, thereby
> allowing iterative access to the file. Otherwise, software
> cannot know that a feature has been fully populated by its
> subfeatures until the end of the file has been reached.
>
>
> =================================================================
>
> REPRESENTING TRANSLATIONS
>
> There are two ways of representing protein translations (e.g. ORFS,
> CDS) in the various implementations of GFF2 and GTF. One way is to
> represent the translation as an interrupted "CDS" region beginning
> with the first base of the first codon and ending at the last base of
> the stop codon. Another is to create a series of exons and to
> indicate the position of the translational start and end on the first
> and last coding exon.
>
> An informal sampling of members of this list (Michele Clamp, Suzi
> Lewis, Richard Durbin) suggests that the latter solution is cleaner
> and more manageable in practice, leading to more consistent annotation
> and to fewer ambiguities. Therefore, I would propose that we
> legislate that translations be represented implicitly by explicit
> translational start and end positions. For this to work properly, the
> parent of the start and end sites must be the mRNA feature and NOT the
> exon.
>
> Under this model, here is a generic gene
>
> gene: a bag of features, including regulatory motifs
> mRNA
> exon
> coding_start
> coding_end
> splice_donor
> splice_acceptor
> 5_utr
> 3_utr
>
> Importantly, the UTRs, coding start and coding end are all children of
> the mRNA. Making them children of the exon (which some will be
> tempted to do!) creates ambiguities in the interpretation of
> alternative splices.
>
> =================================================================
>
> EXAMPLE PROGRAM
>
> I have extended (in an experimental way), the Bio::Tools::GFF module
> to accomodate this new format. Here is a test script and its output
> when run on the above file.
>
> 0 #!/usr/bin/perl -w
> 1 use strict;
> 2 use lib '.';
>
> 3 use Bio::Tools::GFF;
> 4 my $file = 'gff3.txt';
> 5 my $gffio = Bio::Tools::GFF->new(-file=>$file,-gff_version=>3);
> 6 my @f = sort {$a->primary_tag cmp $b->primary_tag} $gffio->features;
> 7 format_features(\@f);
>
> 8 sub format_features {
> 9 my $features = shift;
> 10 my $tabs = shift || 0;
> 11 for my $f (@$features) {
> 12 my $type = $f->primary_tag;
> 13 my $id = $f->unique_id;
> 14 $id ||= '(no id)';# if $id =~ /HASH/;
> 15 my ($start,$end) = ($f->start,$f->end);
> 16 my $hit = $f->can('hstart') ? $f->hunique_id.":".$f->feature2->location->to_FTstring
> 17 : '';
> 18 print "\t"x$tabs,join("\t",$id,$type,$f->location->to_FTstring,$hit),"\n";
> 19 format_features([$f->sub_SeqFeature],$tabs+1);
> 20 }
> 21 }
>
> OUTPUT:
>
> clone00001 clone 1..2679
> contig0001 contig 1..1497228
> gene00001 gene 43733..44677
> mRNA00001 mRNA 43733..44677
> exon00001 exon 43733..43961
> (no id) coding_start 43740
> exon00002 exon 44030..44234
> exon00004 exon 44521..44677
> (no id) coding_end 44677
> cds00001 cds join(43740..43961,44030..44234,44521..44677)
> mRNA00002 mRNA 43733..44677
> exon00001 exon 43733..43961
> (no id) coding_start 43740
> exon00002 exon 44030..44234
> exon00003 exon 44281..44328
> exon00004 exon 44521..44677
> (no id) coding_end 44677
> mRNA03 mRNA 5000..6000
> exon00005 exon 5000..5299
> exon00006 exon 5300..5399
> exon00007 exon 5400..5999
> match0001 match join(1..100,101..500,501..1000,1501..1510) af923:join(1001..1100,1101..1500,1501..2000,2001..2011)
> match0002 match 5001..6000 ua388:1..1000
> (no id) hsp 5001..5500 ua388:1..500
> (no id) hsp 5501..6000 ua388:501..1000
> (no id) repeat 5000..5100
>
>
>
More information about the Bioperl-l
mailing list