<div dir="ltr"><div>Hi all,</div><div><u>The below script parsed BLAST XML output and read a GFF3 file which contains an annotation. If mRNA's ID is in the BLAST output the script will add <i>'Note=Gene description'</i> to mRNA feature.</u></div><div><br></div><div><i>#!/usr/bin/python3<br>import click<br>from Bio.Blast import NCBIXML<br>from BCBio import GFF<br>from pprint import pprint<br><br>class Hits():<br>    def __init__(self, hit_id, hit_def):<br>        self.hit_id = hit_id<br>        self.hit_def = hit_def<br><br>def retrieve_hits_data(blast_XML_file):<br>    with open(blast_XML_file) as bf:<br>        blast_records = NCBIXML.parse(bf)<br>        hits_all = {}<br><br>        for blast_record in blast_records:<br>            query_name = blast_record.query<br>            try:<br>                hit_id = blast_record.alignments[0].hit_id.split('|')[1]<br>                hit_def = blast_record.alignments[0].hit_def.split("OS=")[0]<br>                hits_all[query_name] = Hits(hit_id,hit_def)<br>            except IndexError:<br>                continue<br>    return hits_all<br><br><br>@click.command()<br>@click.option('--gff3', help="Provide GFF3 file", required=True)<br>@click.option('--keep', help="Keep GFF3 file", required=True)<br>@click.option('--reject', help="Reject GFF3 file", required=True)<br>@click.option('--xml', help="Blast XML file", required=True)<br>def run(gff3, keep, reject, xml):<br>    blast_hits = retrieve_hits_data(xml)<br>    # pprint(blast_hits)<br><br>    with open(keep, "w") as out_keep:<br>        with open(reject, "w") as out_reject:<br>            for rec in GFF.parse(gff3):<br>                for feature in rec.features:<br>                    for counter, sub_feature in enumerate(feature.sub_features):<br>                        print(<a href="http://sub_feature.id">sub_feature.id</a>)<br>                        mrna_id = <a href="http://sub_feature.id">sub_feature.id</a><br>                        try:<br>                            blast_hit = blast_hits[mrna_id]<br>                            # feature.qualifiers["Name"] = blast_hit.hit_id<br>                            sub_feature.qualifiers["Note"] = blast_hit.hit_def<br><br>                            print("accepted" + mrna_id + blast_hit.hit_def)<br>                            GFF.write([rec], out_keep)<br>                        except KeyError:<br>                            print("rejected" + mrna_id)<br>                            GFF.write([rec], out_reject)<br>                            continue</i></div><div><i><br></i></div><div><i>if __name__ == '__main__':<br>    run()</i></div><div><br></div><div><u>I failed to extend the above script to only keep mRNAs which have a BLAST hit and the rejected mRNA's features to write to a different file. Below is an example of a GFF3 file and I would like to reject <i>g1</i> and <i>g3</i> because mRNA does not have any BLAST hit. The only <i>g2</i> I would like to keep. What did I miss in the above code?</u><br></div><div><br></div><div><i>##gff-version 3<br>##sequence-region NbV1Ch08 1 129222376<br>NbV1Ch08        AUGUSTUS        <b>gene</b>    7015    29794   0.01    -       .       ID=g1<br>NbV1Ch08        AUGUSTUS        <b>mRNA</b>    7015    29794   0.01    -       .       ID=g1.t1;Parent=g1<br>NbV1Ch08        AUGUSTUS        transcription_end_site  7015    7015    .       -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        three_prime_utr 7015    8531    0.2     -       .       ID=g1.t1.3UTR1;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    7015    8747    .       -       .       ID=g1.t1.exon1;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        stop_codon      8532    8534    .       -       0       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     8532    8747    0.31    -       0       ID=g1.t1.CDS1;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  8748    9191    0.49    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     9192    9342    0.66    -       1       ID=g1.t1.CDS2;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    9192    9342    .       -       .       ID=g1.t1.exon2;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  9343    9915    0.58    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     9916    10006   0.71    -       2       ID=g1.t1.CDS3;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    9916    10006   .       -       .       ID=g1.t1.exon3;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  10007   10101   0.74    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     10102   10201   0.78    -       0       ID=g1.t1.CDS4;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    10102   10201   .       -       .       ID=g1.t1.exon4;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  10202   10712   0.8     -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     10713   11107   0.11    -       2       ID=g1.t1.CDS5;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    10713   11107   .       -       .       ID=g1.t1.exon5;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  11108   11569   0.07    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     11570   12151   0.09    -       2       ID=g1.t1.CDS6;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    11570   12151   .       -       .       ID=g1.t1.exon6;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  12152   12588   0.34    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     12589   12717   0.39    -       2       ID=g1.t1.CDS7;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    12589   12717   .       -       .       ID=g1.t1.exon7;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  12718   12789   0.42    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     12790   13075   0.39    -       0       ID=g1.t1.CDS8;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    12790   13075   .       -       .       ID=g1.t1.exon8;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  13076   14832   0.51    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     14833   15009   0.39    -       0       ID=g1.t1.CDS9;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    14833   15009   .       -       .       ID=g1.t1.exon9;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  15010   15278   0.59    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     15279   15415   0.56    -       2       ID=g1.t1.CDS10;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    15279   15415   .       -       .       ID=g1.t1.exon10;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  15416   15487   0.58    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     15488   15612   0.96    -       1       ID=g1.t1.CDS11;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    15488   15612   .       -       .       ID=g1.t1.exon11;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        intron  15613   15706   0.96    -       .       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        CDS     15707   15957   0.98    -       0       ID=g1.t1.CDS12;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    15707   15958   .       -       .       ID=g1.t1.exon12;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        start_codon     15955   15957   .       -       0       Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        five_prime_utr  15958   15958   0.99    -       .       ID=g1.t1.5UTR1;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        five_prime_utr  27458   28250   0.37    -       .       ID=g1.t1.5UTR2;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    27458   28250   .       -       .       ID=g1.t1.exon13;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        five_prime_utr  29272   29794   0.08    -       .       ID=g1.t1.5UTR3;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        exon    29272   29794   .       -       .       ID=g1.t1.exon14;Parent=g1.t1<br>NbV1Ch08        AUGUSTUS        transcription_start_site        29794   29794   .       -       .       Parent=g1.t1</i></div><div><i>NbV1Ch08        AUGUSTUS        <b>gene</b>    60876   63944   0.03    +       .       ID=g2<br>NbV1Ch08        AUGUSTUS        <b>mRNA</b>    60876   63944   0.03    +       .       ID=g2.t1;<b>Note</b>=B3 domain-containing protein Os03g0120900;Parent=g2<br>NbV1Ch08        AUGUSTUS        transcription_start_site        60876   60876   .       +       .       Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        five_prime_utr  60876   61072   0.19    +       .       ID=g2.t1.5UTR1;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        exon    60876   61072   .       +       .       ID=g2.t1.exon1;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        five_prime_utr  61673   61732   0.37    +       .       ID=g2.t1.5UTR2;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        exon    61673   63449   .       +       .       ID=g2.t1.exon2;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        start_codon     61733   61735   .       +       0       Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        CDS     61733   62974   0.54    +       0       ID=g2.t1.CDS1;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        stop_codon      62972   62974   .       +       0       Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        three_prime_utr 62975   63449   1       +       .       ID=g2.t1.3UTR1;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        three_prime_utr 63565   63944   0.27    +       .       ID=g2.t1.3UTR2;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        exon    63565   63944   .       +       .       ID=g2.t1.exon3;Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        transcription_end_site  63944   63944   .       +       .       Parent=g2.t1<br>NbV1Ch08        AUGUSTUS        <b>gene</b>    64722   65524   0.32    -       .       ID=g3<br>NbV1Ch08        AUGUSTUS        <b>mRNA</b>    64722   65524   0.32    -       .       ID=g3.t1;Parent=g3<br>NbV1Ch08        AUGUSTUS        transcription_end_site  64722   64722   .       -       .       Parent=g3.t1<br>NbV1Ch08        AUGUSTUS        three_prime_utr 64722   64792   0.77    -       .       ID=g3.t1.3UTR1;Parent=g3.t1<br>NbV1Ch08        AUGUSTUS        exon    64722   65524   .       -       .       ID=g3.t1.exon1;Parent=g3.t1<br>NbV1Ch08        AUGUSTUS        stop_codon      64793   64795   .       -       0       Parent=g3.t1<br>NbV1Ch08        AUGUSTUS        CDS     64793   65494   0.44    -       0       ID=g3.t1.CDS1;Parent=g3.t1<br>NbV1Ch08        AUGUSTUS        start_codon     65492   65494   .       -       0       Parent=g3.t1</i></div><div><br></div><div><u>Thank you in advance,</u></div><div><u><br></u></div><div><u>Best wishes,</u></div><div><u><br></u></div><div><u>Michal</u><br></div></div>