2
$\begingroup$

I modified the GFF3 file with the below code:

#!/usr/bin/python3
import click
import gffutils
import gffutils.gffwriter as gffwriter

@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    #db = gffutils.create_db(gff3, dbfn='data/test.db', force=True, keep_order=True, merge_strategy='merge',
    #                        sort_attribute_values=True)

    out_keep = gffwriter.GFFWriter(keep)
    db = gffutils.FeatureDB('data/test.db', keep_order=True)
    for gene in db.features_of_type('gene', order_by='start'):
        print(gene)
        for i in db.children(gene, featuretype='mRNA', order_by='start'):
            print(i)
            i.attributes['Note'] = ["Gene description"]
            print(i)
            out_keep.write_rec(gene)

Unfortunately, the above code saved only this:

#GFF3 file (created by gffutils (v0.9) on 2019-08-28 22:14:13)
NbV1Ch08        AUGUSTUS        gene    7015    29794   0.01    -       .       ID=g1;
NbV1Ch08        AUGUSTUS        gene    60876   63944   0.03    +       .       ID=g2;
NbV1Ch08        AUGUSTUS        gene    64722   65524   0.32    -       .       ID=g3;
NbV1Ch08        AUGUSTUS        gene    66706   67484   0.12    -       .       ID=g4;
NbV1Ch08        AUGUSTUS        gene    115451  122754  0.16    -       .       ID=g5;
NbV1Ch08        AUGUSTUS        gene    147786  155025  0.24    +       .       ID=g6;
NbV1Ch08        AUGUSTUS        gene    155152  161104  0.28    -       .       ID=g7;
NbV1Ch08        AUGUSTUS        gene    237126  238098  0.05    +       .       ID=g8;
NbV1Ch08        AUGUSTUS        gene    265796  297830  0.82    +       .       ID=g9;
NbV1Ch08        AUGUSTUS        gene    265796  297830  0.82    +       .       ID=g9;

Rather than:

NbV1Ch08        AUGUSTUS        gene    7015    29794   0.01    -       .       ID=g1;
NbV1Ch08        AUGUSTUS        mRNA    7015    29794   0.01    -       .       ID=g1.t1;Parent=g1;Note=Gene description;
NbV1Ch08        AUGUSTUS        transcription_end_site  7015    7015    .       -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        three_prime_utr 7015    8531    0.2     -       .       ID=g1.t1.3UTR1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    7015    8747    .       -       .       ID=g1.t1.exon1;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        stop_codon      8532    8534    .       -       0       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     8532    8747    0.31    -       0       ID=g1.t1.CDS1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        intron  8748    9191    0.49    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     9192    9342    0.66    -       1       ID=g1.t1.CDS2;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    9192    9342    .       -       .       ID=g1.t1.exon2;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  9343    9915    0.58    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     9916    10006   0.71    -       2       ID=g1.t1.CDS3;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    9916    10006   .       -       .       ID=g1.t1.exon3;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  10007   10101   0.74    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     10102   10201   0.78    -       0       ID=g1.t1.CDS4;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    10102   10201   .       -       .       ID=g1.t1.exon4;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  10202   10712   0.8     -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     10713   11107   0.11    -       2       ID=g1.t1.CDS5;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    10713   11107   .       -       .       ID=g1.t1.exon5;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  11108   11569   0.07    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     11570   12151   0.09    -       2       ID=g1.t1.CDS6;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    11570   12151   .       -       .       ID=g1.t1.exon6;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  12152   12588   0.34    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     12589   12717   0.39    -       2       ID=g1.t1.CDS7;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    12589   12717   .       -       .       ID=g1.t1.exon7;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  12718   12789   0.42    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     12790   13075   0.39    -       0       ID=g1.t1.CDS8;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    12790   13075   .       -       .       ID=g1.t1.exon8;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  13076   14832   0.51    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     14833   15009   0.39    -       0       ID=g1.t1.CDS9;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    14833   15009   .       -       .       ID=g1.t1.exon9;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15010   15278   0.59    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15279   15415   0.56    -       2       ID=g1.t1.CDS10;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15279   15415   .       -       .       ID=g1.t1.exon10;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15416   15487   0.58    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15488   15612   0.96    -       1       ID=g1.t1.CDS11;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15488   15612   .       -       .       ID=g1.t1.exon11;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        intron  15613   15706   0.96    -       .       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        CDS     15707   15957   0.98    -       0       ID=g1.t1.CDS12;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    15707   15958   .       -       .       ID=g1.t1.exon12;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        start_codon     15955   15957   .       -       0       Parent=g1.t1;
NbV1Ch08        AUGUSTUS        five_prime_utr  15958   15958   0.99    -       .       ID=g1.t1.5UTR1;Parent=g1.t1
NbV1Ch08        AUGUSTUS        five_prime_utr  27458   28250   0.37    -       .       ID=g1.t1.5UTR2;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    27458   28250   .       -       .       ID=g1.t1.exon13;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        five_prime_utr  29272   29794   0.08    -       .       ID=g1.t1.5UTR3;Parent=g1.t1
NbV1Ch08        AUGUSTUS        exon    29272   29794   .       -       .       ID=g1.t1.exon14;Parent=g1.t1;
NbV1Ch08        AUGUSTUS        transcription_start_site        29794   29794   .       -       .       Parent=g1.t1;

What did I miss?

Thank you in advance,

Update This is base on @daniel-standage code.

import click
import tag
def add_description(instream):
    for feature in instream:
        for mrna in tag.select.features(feature, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')

        yield feature


@click.command()
@click.option('--gff3', help="Provide GFF3 file", required=True)
@click.option('--keep', help="Keep GFF3 file", required=True)
@click.option('--reject', help="Reject GFF3 file", required=True)
def run(gff3, keep, reject):
    reader = tag.GFF3Reader(infilename=gff3)
    desc_stream = add_description(reader)
    writer = tag.GFF3Writer(desc_stream)
    writer.write()


if __name__ == '__main__':
    run()

which caused this error:

mod_braker_gff3.py --gff3 data/NbV1Ch08-augustus.hints_utr.gff3 --keep data/keep.gff3 --reject data/reject.gff3
Traceback (most recent call last):
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 26, in <module>
    run()
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 764, in __call__
    return self.main(*args, **kwargs)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 717, in main
    rv = self.invoke(ctx)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 956, in invoke
    return ctx.invoke(self.callback, **ctx.params)
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/click/core.py", line 555, in invoke
    return callback(*args, **kwargs)
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 22, in run
    writer.write()
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/writer.py", line 57, in write
    for entry in self._instream:
  File "/Users/lorencm/projects/bioinf-scripts/mod_braker_gff3.py", line 8, in add_description
    for mrna in tag.select.features(feature, type='mRNA', traverse=True):
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/select.py", line 24, in features
    for feature in entry_type_filter(entrystream, tag.Feature):
  File "/anaconda/envs/bioinf-scripts/lib/python3.7/site-packages/tag/select.py", line 99, in entry_type_filter
    for entry in entrystream:
TypeError: 'Directive' object is not iterable

Process finished with exit code 1
| improve this question | |
$\endgroup$
  • $\begingroup$ As haci, said you're filtering for genes. When you have this sort of issue, the first step should always be to print everything. For example, if you tried print db.children(gene, featuretype='mRNA', order_by='start'), you'd presumably see that it is empty. $\endgroup$ – terdon Aug 28 '19 at 13:01
  • 2
    $\begingroup$ If you cross post your question please leave a note, so that other can see what is already answered! $\endgroup$ – finswimmer Aug 29 '19 at 4:33
2
$\begingroup$

I guess you want to take the features for a gene?

As @haci and @terdon pointed out, it's a problem with the selected. Below I included an example, with mouse_extra_comma.gff3 obtained from example file provided by the package gffutils.

#gff.py
#!/usr/bin/python3

###
def addNote(x):
    x.attributes['Note'] = ["Gene description"]
    return(x)

###

import click
import gffutils
import gffutils.gffwriter as gffwriter

fn = gffutils.example_filename('mouse_extra_comma.gff3')
db = gffutils.create_db(fn, dbfn='test.db', force=True, keep_order=True,
    merge_strategy='merge', sort_attribute_values=True)

db = gffutils.FeatureDB('./test.db', keep_order=True)

## here's to test so you can see how it works
gene = next(db.features_of_type('gene', order_by='start')) 
items =[addNote(i) for i in db.children(gene, featuretype=['exon','CDS'],order_by='start')]

items[0:5]
#[<Feature mRNA (chr17:6797760-6818159[+]) at 0x103659c50>,
# <Feature exon (chr17:6797760-6797769[+]) at 0x1037cee10>,
# <Feature five_prime_UTR (chr17:6797760-6797769[+]) at 0x103660898>,
# <Feature exon (chr17:6806513-6806553[+]) at 0x103659f98>,
# <Feature five_prime_UTR (chr17:6806513-6806526[+]) at 0x1036595f8>]
#####
## here's following your... code...
for gene in db.features_of_type('gene', order_by='start'):
    print(gene)
    for i in db.children(gene,order_by='start'):
        print(addNote(i))

Once you define the feature you want, should be ok

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you but I do not know how to do it? $\endgroup$ – user977828 Aug 29 '19 at 0:02
  • $\begingroup$ Could you elaborate on which step are you stuck? Would be glad to guide you along $\endgroup$ – StupidWolf Aug 29 '19 at 7:58
  • $\begingroup$ Where could I add something like `mrna.attributes['Note'] = "Gene description"? $\endgroup$ – user977828 Aug 29 '19 at 14:12
  • $\begingroup$ You would write a function. I have edited my code above, although i am still a bit confused about what you need. If I get you correct, you would like to: 1. For every gene, print the gene first 2. Then print all features (mRNA, exon, CDS) associated with the gene, and attached a gene description note 3. Write to this a file Is this correct? If so, could you please formulate the question in this way. $\endgroup$ – StupidWolf Aug 29 '19 at 19:18
  • $\begingroup$ Thank you for the code. I need only to append to the original GFF3 file in the mRNA attribute Note=Gene description. For example from NbV1Ch08 AUGUSTUS mRNA 22418839 22429909 0.17 + . ID=g1168.t1;Parent=g1168; to NbV1Ch08 AUGUSTUS mRNA 22418839 22429909 0.17 + . ID=g1168.t1;Parent=g1168;Note=Gene description;. Next I would like to save the changes. $\endgroup$ – user977828 Aug 29 '19 at 23:05
1
$\begingroup$

You might want to check this line out:

for gene in db.features_of_type('gene', order_by='start'):

You are specifically filtering for genes hence other types of annotations are filtered out.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you but I do not know how to do it? $\endgroup$ – user977828 Aug 29 '19 at 0:03
  • $\begingroup$ Please see this about the features_of_type method: pythonhosted.org/gffutils/autodocs/… . It basically tells to provide a string or tuple of feature(s) that you would like to retrieve, in your case all the features in the original file. The easiest would be take the unique values in the "feature" column of your input and feed it into the "features_of_type" method. $\endgroup$ – haci Aug 29 '19 at 0:15
  • $\begingroup$ Where could I add something like `mrna.attributes['Note'] = "Gene description"? $\endgroup$ – user977828 Aug 29 '19 at 14:15
1
$\begingroup$

If you're willing to give another GFF3 processing library a shot, I'd recommend checking out the tag library. The following code should add the note attribute to mRNA records as you're attempting to do in you're original code, but will also print out all other records.

import sys
import tag

def add_description(instream):
    for record in instream:
        if not isinstance(record, tag.Feature):
            yield record
            continue
        for mrna in tag.select.features(record, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')
        yield feature

reader = tag.GFF3Reader(instream=sys.stdin)
desc_stream = add_description(reader)
writer = tag.GFF3Writer(desc_stream)
writer.retainids = True
writer.write()

I didn't use click or argparse to put any CLI on the script, so unless you want to add those just use STDIN and STDOUT for I/O.

$ python3 myscript.py < mydata.gff3 > mydatawithdescription.gff3

P.S. The add_description can be simplified if you don't care about retaining comments, ##sequence-region pragmas, and so on.

def add_description(instream):
    for feature in tag.select.features(instream):
        for mrna in tag.select.features(record, type='mRNA', traverse=True):
            mrna.add_attribute('Note', 'Gene description')
        yield feature
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your solution. Unfortunately, I ran into TypeError: 'Directive' object is not iterable as shown in my updated question. What did I miss? $\endgroup$ – user977828 Aug 29 '19 at 0:01
  • $\begingroup$ @user977828 I wrote the initial code with the (bad) assumption that the input data didn't include any comments, pragmas, directives, etc. I've updated my answer to show how to handle input that includes these types of records, or below that is a simpler solution that ignores them. $\endgroup$ – Daniel Standage Aug 29 '19 at 18:29
  • $\begingroup$ Thank you it works. $\endgroup$ – user977828 Aug 29 '19 at 23:14
  • $\begingroup$ How is it possible to avoid to change mRNA's ID? $\endgroup$ – user977828 Aug 30 '19 at 12:44
  • $\begingroup$ At the moment, there is no way to change tag's default behavior: it rewrites all ID attributes, since according to the spec these attributes are only guaranteed to be consistent and unique within a single file, not across files. $\endgroup$ – Daniel Standage Aug 30 '19 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.