GFF3 and assembly handling, TAIR12 - wdlingit/cop GitHub Wiki

This document describes my process of handling TAIR12 GFF3 file at date 20260910. Refer the environment setting if needed.

$ date
Thu Sep 10 11:32:23 CST 2026

$ ls -l
total 80523
-rwxrwxrwx 1 ubuntu ubuntu 16465862 Sep 10 14:29 Araport11_GFF3_genes_transposons.20250813.gff.gz
-rwxrwxrwx 1 ubuntu ubuntu 39229219 Sep 10 14:29 Col-CC_v2_genome.fasta.gz
-rwxrwxrwx 1 ubuntu ubuntu 37482399 Jun 17  2023 GCF_000001735.4_TAIR10.1_genomic.fna.gz
-rwxrwxrwx 1 ubuntu ubuntu   156517 Sep 10 14:32 TAIR10_chrC.fas
-rwxrwxrwx 1 ubuntu ubuntu   371653 Sep 10 14:32 TAIR10_chrM.fas
-rwxrwxrwx 1 ubuntu ubuntu 11093467 Sep 10 14:29 TAIR12_1Feb26.gff3.gz

File descriptions:

Col-CC_v2_genome.fasta.gz: TAIR12 genome assembly, contains Chr1~5
TAIR12_1Feb26.gff3.gz: TAIR12 genome annotation
Araport11_GFF3_genes_transposons.20250813.gff.gz: Newest Araport11 genome annotation in the TAIR dabase
GCF_000001735.4_TAIR10.1_genomic.fna.gz: One possible source for ChrC and ChrM, from NCBI TAIR10.1
TAIR10_chrC.fas and TAIR10_chrM.fas: Another possible source for ChrC and ChrM, from TAIR

Quote from README_whole_chromosomes.txt in TAIR.

--- Organellar chromosomes (used by TAIR12, unchanged) ---

TAIR12 did not re-assemble or re-curate the chloroplast or mitochondrial genomes.
They remain the TAIR10 sequences, in TAIR10/:

  TAIR10/TAIR10_chrC.fas   (chloroplast)
  TAIR10/TAIR10_chrM.fas   (mitochondria)

Quote from README of TAIR12.

Mitochondrial and chloroplast genomes were not updated or re-curated in TAIR12-- their gene models remain identical to those in Araport11.

Initial parse of GFF3

Some preprocessing of adding ChrC and ChrM.

# extract ChrC ChrM part from Araport11 genome annotation
$ cat Araport11_GFF3_genes_transposons.20250813.gff | perl -ne 'print if /^ChrC/ || /^ChrM/' > Araport11_GFF3.chrCchrM.gff

# put Chr1~5/C/M genome annotation together into file TAIR12_1Feb26.ChrC_ChrM.gff3
$ cp TAIR12_1Feb26.gff3 TAIR12_1Feb26.ChrC_ChrM.gff3
$ cat Araport11_GFF3.chrCchrM.gff >> TAIR12_1Feb26.ChrC_ChrM.gff3

It was found that some CDS records are sharing the same IDs. The following command was applied to removing those ID's.

$ cat TAIR12_1Feb26.ChrC_ChrM.gff3 | perl -ne 'if(/^#/){ print }else{ chomp; @t=split(/\t/); if($t[2] eq "CDS"){ $t[8]=~s/ID=.+?;// } print join("\t",@t)."\n"; }' > TAIR12_1Feb26.ChrC_ChrM.fixed.gff

It was also found that two tRNA's got no exon records, and they were manually added.

$ diff TAIR12_1Feb26.ChrC_ChrM.fixed.gff TAIR12_1Feb26.ChrC_ChrM.fixed1.gff
797357a797358
> ChrM  Araport11       exon    149673  149743  .       +       .       Parent=ATMG00030.1;
797456a797458
> ChrM  Araport11       exon    250307  250380  .       +       .       Parent=ATMG00460.1;

Applied misc.GffTree to check and parse the GFF3 file.

$ java misc.GffTree -I TAIR12_1Feb26.ChrC_ChrM.fixed1.gff
ID attribute (-idAttr) not assigned, using default: ID
parent attribute list (-parentAttr) not assigned, using default: [Parent, Derives_from]
program: GffTree
input GFF file (-I): TAIR12_1Feb26.ChrC_ChrM.fixed1.gff
ID attribute (-idAttr): ID
parent attribute list (-parentAttr): [Parent, Derives_from]

# feature structure inside the GFF3 file
wdlin@comp01:/RAID2/R418/20260831_Grillet/DB/tair12_ex$ cat TAIR12_1Feb26.ChrC_ChrM.fixed1.gff.features
GffRoot
  repeat_region*
    repeat_region*
    mobile_element*
  gene*
    mRNA*
      exon
      CDS
      misc_feature*
      protein*
    ncRNA*
      exon
    tRNA*
      exon
    transcript*
      exon
      misc_feature*
    miRNA_primary_transcript*
      exon
    gene*
      ncRNA*
        exon
    rRNA*
      exon
    precursor_RNA*
    misc_feature*
    misc_RNA*
      exon
    uORF*
    exon
  mobile_element*
    mobile_element*
    repeat_region*
  pseudogene*
    transcript*
      exon
      misc_feature*
    misc_RNA*
      exon
    misc_feature*
    pseudogenic_tRNA*
      pseudogenic_exon*

Second-level feature extraction

Since there are gene records that might be mRNA, rRNA, ncRNA... etc (at second level of the feature tree). It is desirable to classify gene categories following the info inside the GFF3 file. Firstly, second-level features were listed in second.features.

$ cat second.features
mRNA
ncRNA
tRNA
transcript
miRNA_primary_transcript
rRNA
misc_RNA
transcript
misc_RNA
pseudogenic_tRNA

However, some features like misc_RNA and transcript may not be desirable name of categories because they belong to first-level features gene and pseudogene at the same time. So for misc_RNA and transcript records actually under pseudogene, pseudogene should be the appropriate category name. Accordingly, the category assienment process should take features of parents into consideration.

$ cat TAIR12_1Feb26.ChrC_ChrM.fixed1.gff | perl -ne '
    if($.==1){
        open(FILE,"<second.features");
        while($line=<FILE>){
            chomp $line;
            $hash{$line}=1
        }
        close FILE;
    }
    chomp;
    @t=split(/\t/);
    $t[8].=";";
    if($t[8]=~/ID=(.+?);/){
        $feature0{$1}=$t[2];
    }
    if(exists $hash{$t[2]}){
        $t[8]=~/Parent=(.+?);/ || $t[8]=~/Derives_from=(.+?);/;
        $feature{$1}{$t[2]}=1
    }
    if(eof STDIN){
        print "geneID\tfeature\n";
        for $k (sort keys %feature){
            $feature{$k}{$feature0{$k}}=1 if exists $feature0{$k};
            @features=sort keys %{$feature{$k}};
            print "$k\t@features\n";
        }
    }
' | perl -ne '
    chomp;
    @t=split;
    $id=shift @t;
    if($t[0] eq "gene"){ shift @t }
    if(($t[0] eq "pseudogene") && ($t[1] eq "transcript")){ pop @t }
    if(($t[0] eq "pseudogene") && @t>1){ shift @t }
    if(($t[0] eq "misc_RNA") && @t>1){ shift @t }
    if(($t[0] eq "miRNA_primary_transcript") && @t>1){ pop @t }
    print "$id\t@t\n";
' > gene.features

Explanation of the above command:

The first perl oneliner was to collect feature of every parent, as well as their children's feature.
The second perl oneliner was to do appropriate postprocessing. For example, if a gene came with gene and ncRNA, we should take ncRNA but not gene.

Interval extraction

Applied misc.ModelCGFF for interval coordinate extraction. Note that the -GRE triple for miRNA_primary_transcript is -GRE gene miRNA_primary_transcript miRNA_primary_transcript. The -GRE triple for rRNA is -GRE gene rRNA rRNA.

$ java misc.ModelCGFF -GFF3 TAIR12_1Feb26.ChrC_ChrM.fixed1.gff -GRE gene mRNA:ncRNA:tRNA:transcript:misc_RNA exon:CDS -GRE gene miRNA_primary_transcript miRNA_primary_transcript -GRE gene rRNA rRNA -GRE pseudogene transcript:misc_RNA:pseudogenic_tRNA exon:pseudogenic_exon -O TAIR12.strand -IP true
ID attribute (-idAttr) not assigned, using default: ID
parent attribute list (-parentAttr) not assigned, using default: [Parent, Derives_from]
program: ModelCGFF
GFF3 filename (-GFF3): TAIR12_1Feb26.ChrC_ChrM.fixed1.gff
ID attribute (-idAttr): ID
parent attribute list (-parentAttr)): [Parent, Derives_from]
gene-rna-exon feature triples (-GRE): [[gene, mRNA:ncRNA:tRNA:transcript:misc_RNA, exon:CDS], [gene, miRNA_primary_transcript, miRNA_primary_transcript], [gene, rRNA, rRNA], [pseudogene, transcript:misc_RNA:pseudogenic_tRNA, exon:pseudogenic_exon]]
intron preserving merged model (-IP)): true
representative list (-rep): null
output prefix (-O): TAIR12.strand

wdlin@comp01:/RAID2/R418/20260831_Grillet/DB/tair12_ex$ head TAIR12.strand.cgff
>AT1G01010      Chr1    7395    9666    +
7395    7677
7760    8040
8250    8369
8470    8859
8938    9090
9203    9666
>AT1G01020      Chr1    10572   12478   -
10572   10834
10922   10997
wdlin@comp01:/RAID2/R418/20260831_Grillet/DB/tair12_ex$ head TAIR12.strand.model
>AT1G01010.1    Chr1    7395    9666    +       AT1G01010
7395    7677
7760    8040
8250    8369
8470    8859
8938    9090
9203    9666
>AT1G01020.3    Chr1    10572   12458   -       AT1G01020
10572   10834
10922   10997

Gene info extraction

File features lists interested features and their level in the feature tree.

$ cat features
1       gene
2       mRNA
2       ncRNA
2       tRNA
2       transcript
2       miRNA_primary_transcript
2       rRNA
2       misc_RNA
1       pseudogene
2       transcript
2       misc_RNA
2       pseudogenic_tRNA

Numbers of attributes associated to level-1 features and level-2 features were observed.

$ cat TAIR12_1Feb26.ChrC_ChrM.fixed1.gff | perl -ne 'if($.==1){ open(FILE,"<features"); while($line=<FILE>){ chomp $line; @s=split(/\s+/,$line); $feature{$s[1]}=$s[0]; } close FILE; } chomp; @t=split(/\t/); if($feature{$t[2]}==1){ $t[8].=";" if $t[8]!~/;$/; $tmpstr=$t[8]; while($tmpstr=~/^([^,;=]+?)=(.+)$/){ $key=$1; $hash{$key}++; $tmpstr=$2; if($tmpstr=~/(.+?);\s*([\S^=]+?=.+)/){ $val=$1; $tmpstr=$2 }else{ ($val)=$tmpstr=~/(.+);/; $tmpstr=""; } } } if(eof){ for $k (sort keys %hash){ print "$k\t$hash{$k}\n" } }'
Alias   24
Dbxref  40369
Derives_from    68
ID      40505
Name    40504
Note    5022
computational_description       125
copy_num_id     36
coverage        179
curator_summary 123
description     13375
evidence_code   2
exception       2
extra_copy_number       179
family_EDS      3994
full_name       64
gbkey   6
gene    16014
gene_synonym    7281
id2     6
locus   85
locus_biotype   40236
locus_tag       40236
locus_type      267
other   117
part    2
product 889
pseudo  3
sequence_id     179
symbol  131

$ cat TAIR12_1Feb26.ChrC_ChrM.fixed1.gff | perl -ne 'if($.==1){ open(FILE,"<features"); while($line=<FILE>){ chomp $line; @s=split(/\s+/,$line); $feature{$s[1]}=$s[0]; } close FILE; } chomp; @t=split(/\t/); if($feature{$t[2]}==2){ $t[8].=";" if $t[8]!~/;$/; $tmpstr=$t[8]; while($tmpstr=~/^([^,;=]+?)=(.+)$/){ $key=$1; $hash{$key}++; $tmpstr=$2; if($tmpstr=~/(.+?);\s*([\S^=]+?=.+)/){ $val=$1; $tmpstr=$2 }else{ ($val)=$tmpstr=~/(.+);/; $tmpstr=""; } } } if(eof){ for $k (sort keys %hash){ print "$k\t$hash{$k}\n" } }'
Dbxref  2708
ID      54027
Name    54064
Note    4021
Parent  54027
computational_description       190
conf_class      1
conf_rating     88
confidence      5664
contains        231
cpc2_call       5485
cpc2_score      5485
curator_summary 188
extra_copy_number       179
full_name       64
locus_biotype   2645
locus_tag       2645
locus_type      33
longread_counts 5485
ncRNA_class     7279
num_assemblies  5485
orf_len 5485
pre_miRNA       102
product 47097
rfam_acc        164
rfam_id 164
riboseq 476
symbol  190
tau     5485

For this time, it was decided to extract attributes from level-1 features.

$ cat TAIR12_1Feb26.ChrC_ChrM.fixed1.gff | perl -ne '
    if($.==1){ 
        open(FILE,"<features"); 
        while($line=<FILE>){ 
            chomp $line; 
            @s=split(/\s+/,$line); 
            $feature{$s[1]}=$s[0]; 
        } 
        close FILE; 
    } 
    chomp; @t=split(/\t/); 
    if($feature{$t[2]}==1){ 
        $t[8].=";" if $t[8]!~/;$/; 
        ($id)=$t[8]=~/ID=(.+?);/; 
        $tmpstr=$t[8]; 
        while($tmpstr=~/^([^,;=]+?)=(.+)$/){ 
            $key=$1; 
            $tmpstr=$2; 
            if($tmpstr=~/(.+?);\s*([\S^=]+?=.+)/){ 
                $val=$1; 
                $tmpstr=$2 
            }else{ 
                ($val)=$tmpstr=~/(.+);/; 
            } 
            $hash{$id}{$key}=$val; 
        } 
    } 
    if(eof STDIN){ 
        print "#transcript\tNote\tdescription\n"; 
        for $id (sort keys %hash){ 
            print "$id"; 
            for $k ("Note","description"){ 
                if(exists $hash{$id}{$k}){ 
                    print "\t$hash{$id}{$k}" 
                }else{ 
                    print "\t" 
                } 
            } 
            print "\n" 
        } 
    }
' > TAIR12.description.txt

The following process is the my current best practice to fix unprocessable text. Firstly, visually confirm the following code decodes those encoded text, as there might be some unprocessable spots.

$ cat TAIR12.description.txt | perl -MHTML::Entities -Xne '
    while(/(\S+\s+)/g){ 
        $x=$1; 
        $x=~s/%([0-9A-Fa-f]{2})/chr(hex($1))/eg; 
        print decode_entities($x) 
    }
' | perl -ne '
    print if /[^[:ascii:]]/
' | less

If there is any line with unprocessable characters, it is needed to manually correct them. For this time, there is no unprocessable characters inside the file. So we don't have to do the manual correction. However, some info was actually encoded and we have to decode them.

$ cat TAIR12.description.txt | perl -MHTML::Entities -Xne '
    while(/(\S+\s+)/g){ 
        $x=$1; 
        $x=~s/%([0-9A-Fa-f]{2})/chr(hex($1))/eg; 
        print decode_entities($x) 
    }
' > TAIR12.description.decoded.txt

Check: TAIR12 Chr1~5 + NCBI TAIR10.1 ChrC/ChrM

Manually edit Col-CC_v2_genome.fasta and add ChrC and ChrM from GCF_000001735.4_TAIR10.1_genomic.fna and save them in TAIR12.genome.fa.

$ fixFasta.pl TAIR12.genome.fa TAIR12.genome.fixed.fa

Extract CDS coordinates for sequence generation.

$ java misc.CanonicalGFF -I TAIR12_1Feb26.ChrC_ChrM.fixed1.gff -O TAIR12.coding.cgff -GE mRNA CDS
ID attribute (-idAttr) not assigned, using default: ID
parent attribute list (-parentAttr) not assigned, using default: [Parent, Derives_from]
program: CanonicalGFF
input GFF file (-I): TAIR12_1Feb26.ChrC_ChrM.fixed1.gff
output filename (-O): TAIR12.coding.cgff
gene-exon feature pairs (-GE): [[mRNA, CDS]]
ID attribute (-idAttr): ID
parent attribute list (-parentAttr)): [Parent, Derives_from]
intron preserving (-IP)): false

Generate CDS sequences and translate them.

$ SeqGen.pl TAIR12.coding.fasta TAIR12.genome.fixed.fa TAIR12.coding.cgff

$ translate.pl TAIR12.coding.fasta TAIR12.coding.faa

Check translated sequences by checking Stop signs inside protein sequences.

$ cat TAIR12.coding.faa | perl -ne 'chomp; if(/^>(\S+)/){ $id=$1 }else{ print "$id\n" if /\*./ }'
ATMG01275.1
ATMG00580.1
ATMG00650.1

Check: TAIR12 Chr1~5 + TAIR TAIR10 ChrC/ChrM

Manually edit Col-CC_v2_genome.fasta and add TAIR10 ChrC and TAIR10 ChrM from the TAIR database and save them in TAIR12.test.fa.

$ fixFasta.pl TAIR12.test.fa TAIR12.test.fixed.fa

We applied the same extracted CDS coordinates for generating CDS sequences and translate them.

$ SeqGen.pl TAIR12.test.coding.fasta TAIR12.test.fixed.fa TAIR12.coding.cgff

$ translate.pl TAIR12.test.coding.fasta TAIR12.test.coding.faa