Sequence annotation

Gene models are read from GFF3 (or GTF) into a SequenceAnnotation: a directed acyclic graph of SequenceInterval objects (genes, mRNAs, exons, CDSs, …) linked by their Parent attributes.

from picea import SequenceAnnotation, SequenceInterval

Single intervals

One GFF3 line corresponds to one interval. The first eight columns become attributes, and so does every key in the ninth column. Attribute keys are lowercased (except ID) and attribute values are always lists.

interval = SequenceInterval.from_gff_line("ctg123\t.\tgene\t1000\t9000\t.\t+\t.\tID=gene00001;Name=EDEN")
interval

<SequenceInterval type=gene ID=gene00001 loc=ctg123..1000..9000..+ at 0x7f0a789fc050>

interval.seqid, interval.start, interval.end, interval.strand, interval.name

('ctg123', 1000, 9000, '+', ['EDEN'])

interval.to_gff_line()

'ctg123\t.\tgene\t1000\t9000\t.\t+\t.\tID=gene00001;Name=EDEN'

Gene models

This file contains a single gene from the Medicago truncatula genome annotation.

annotation = SequenceAnnotation.from_gff(filename="data/genemodel.gff3")
len(annotation)

21

Intervals are accessed by ID, and can be grouped or filtered with any function of an interval.

by_type = annotation.groupby(lambda interval: interval.interval_type)
{interval_type: len(intervals) for interval_type, intervals in by_type.items()}

{'gene': 1,
 'mRNA': 1,
 'five_prime_UTR': 1,
 'exon': 8,
 'CDS': 8,
 'repeat_region': 1,
 'three_prime_UTR': 1}

long_exons = annotation.filter(lambda interval: interval.interval_type == "exon" and interval.end - interval.start > 500)
[(exon.ID, exon.end - exon.start) for exon in long_exons]

[('exon:MtrunA17Chr1g0184451.8', 574)]

children and parents are transitive: they contain all descendants and all ancestors of an interval. They return a new SequenceAnnotation, so groupby and filter work on them as well.

gene = annotation["gene:MtrunA17Chr1g0184451"]
for child in gene.children:
    print(f"{child.interval_type:16}{child.start:>10}{child.end:>10}  {child.ID}")

mRNA              35039693  35045433  mRNA:MtrunA17Chr1g0184451
five_prime_UTR    35039693  35039922  five_prime_UTR:MtrunA17Chr1g0184451.0
exon              35039693  35039934  exon:MtrunA17Chr1g0184451.1
CDS               35039923  35039934  CDS:MtrunA17Chr1g0184451.1
exon              35040034  35040142  exon:MtrunA17Chr1g0184451.2
CDS               35040034  35040142  CDS:MtrunA17Chr1g0184451.2
exon              35041155  35041277  exon:MtrunA17Chr1g0184451.3
CDS               35041155  35041277  CDS:MtrunA17Chr1g0184451.3
exon              35042183  35042280  exon:MtrunA17Chr1g0184451.4
CDS               35042183  35042280  CDS:MtrunA17Chr1g0184451.4
exon              35042442  35042568  exon:MtrunA17Chr1g0184451.5
CDS               35042442  35042568  CDS:MtrunA17Chr1g0184451.5
exon              35043200  35043285  exon:MtrunA17Chr1g0184451.6
CDS               35043200  35043285  CDS:MtrunA17Chr1g0184451.6
exon              35044019  35044144  exon:MtrunA17Chr1g0184451.7
CDS               35044019  35044144  CDS:MtrunA17Chr1g0184451.7
CDS               35044859  35045017  CDS:MtrunA17Chr1g0184451.8
exon              35044859  35045433  exon:MtrunA17Chr1g0184451.8
three_prime_UTR   35045018  35045433  three_prime_UTR:MtrunA17Chr1g0184451.16

cds = by_type["CDS"].elements[0]
[(parent.interval_type, parent.ID) for parent in cds.parents]

[('mRNA', 'mRNA:MtrunA17Chr1g0184451'), ('gene', 'gene:MtrunA17Chr1g0184451')]

gene.children.groupby(lambda interval: interval.interval_type)["exon"].elements

[<SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.1 loc=MtrunA17Chr1..35039693..35039934..+ at 0x7f0a78808c30>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.2 loc=MtrunA17Chr1..35040034..35040142..+ at 0x7f0a78808e90>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.3 loc=MtrunA17Chr1..35041155..35041277..+ at 0x7f0a78809220>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.4 loc=MtrunA17Chr1..35042183..35042280..+ at 0x7f0a78809480>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.5 loc=MtrunA17Chr1..35042442..35042568..+ at 0x7f0a788096e0>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.6 loc=MtrunA17Chr1..35043200..35043285..+ at 0x7f0a78809940>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.7 loc=MtrunA17Chr1..35044019..35044144..+ at 0x7f0a78809ba0>,
 <SequenceInterval type=exon ID=exon:MtrunA17Chr1g0184451.8 loc=MtrunA17Chr1..35044859..35045433..+ at 0x7f0a78809f30>]

Output

Intervals and annotations can be written to GFF3, GTF, and json.

gene.to_gff_line()

'MtrunA17Chr1\tEuGene\tgene\t35039693\t35045433\t.\t+\t.\tID=gene:MtrunA17Chr1g0184451;Name=MtrunA17Chr1g0184451;locus_tag=MtrunA17_Chr1g0184451'

print(gene.to_json(indent=2))

{
  "ID": "gene:MtrunA17Chr1g0184451",
  "seqid": "MtrunA17Chr1",
  "source": "EuGene",
  "interval_type": "gene",
  "start": 35039693,
  "end": 35045433,
  "score": ".",
  "strand": "+",
  "phase": ".",
  "attributes": {
    "name": [
      "MtrunA17Chr1g0184451"
    ],
    "locus_tag": [
      "MtrunA17_Chr1g0184451"
    ]
  }
}

GTF files have no interval IDs or Parent attributes: intervals are linked by their gene_id and transcript_id attributes instead. Reading GTF recreates the gene model.

gtf = annotation.to_gtf()
print("\n".join(gtf.split("\n")[:3]))

MtrunA17Chr1	EuGene	gene	35039693	35045433	.	+	.	name "MtrunA17Chr1g0184451"; locus_tag "MtrunA17_Chr1g0184451"; gene_id "gene:MtrunA17Chr1g0184451";
MtrunA17Chr1	EuGene	transcript	35039693	35045433	.	+	.	name "MtrunA17Chr1g0184451"; ontology_term "GO:0005739"; locus_tag "MtrunA17_Chr1g0184451"; product "Putative flagellum site-determining protein YlxH/ Fe-S cluster assembling factor NBP35"; transcript_id "mRNA:MtrunA17Chr1g0184451"; gene_id "gene:MtrunA17Chr1g0184451";
MtrunA17Chr1	EuGene	five_prime_UTR	35039693	35039922	.	+	.	est_cons "100.0"; est_incons "0.0"; transcript_id "mRNA:MtrunA17Chr1g0184451"; gene_id "gene:MtrunA17Chr1g0184451";

from_gtf = SequenceAnnotation.from_gtf(string=gtf)
[(interval.interval_type, interval.ID) for interval in from_gtf["gene:MtrunA17Chr1g0184451"].children][:4]

[('mRNA', 'mRNA:MtrunA17Chr1g0184451'),
 ('five_prime_UTR', 'mRNA:MtrunA17Chr1g0184451.five_prime_UTR_0'),
 ('exon', 'mRNA:MtrunA17Chr1g0184451.exon_0'),
 ('CDS', 'mRNA:MtrunA17Chr1g0184451.CDS_0')]