<p>In your case, I would definitely suggest following @Emily_Ensembl's advice and using the Arabidopsis GTF from Ensembl. But for future reference, if an Ensembl GTF wasn't available, you could build something like this using the <code>gtf</code> class from <code>cgat</code></p>

<p>The cgat module and dependancies can be installed by following the instructions <a href="https://github.com/cgat-developers/cgat-apps" rel="nofollow noreferrer">here</a>. Specifically, you can download and run their installation script. This will set up a dedicated conda environment to run the versions of dependancies needed. Further documentation can be found <a href="https://www.cgat.org/downloads/public/cgat/documentation/" rel="nofollow noreferrer">here</a>.</p>

<pre><code># download installation script:
curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/master/install.sh

# install the development version (recommended, no production version yet):
bash install.sh --devel
</code></pre>

<p>This requires python3 and anaconda to install. You may need to run this script in the conda environment:</p>

<pre><code># show available environments
conda info --envs
# activate by path
source activate /home/user/local/bin/conda-install/envs/cgat-a
</code></pre>

<p>Then run the follow as a python3 script. Copy the contents in a file such as <code>convert_gtf.py</code>. Then run it in the terminal on your gtf file:</p>

<pre><code>python convert_gtf.py genes.gtf &gt; genes_new.gtf
</code></pre>

<pre class="lang-python prettyprint-override"><code>#python 3.6.4

import sys

from cgat import GTF
from cgatcore import iotools

for gene in GTF.flat_gene_iterator(
     GTF.iterator(iotools.open_file(sys.argv[1])), strict = False):
    gene_start = min (e.start for e in gene)
    gene_end = max(e.end for e in gene)

    if any(e.feature == "CDS" for e in gene):
        gene_type = "protein_coding"
    else:
        gene_type = "non_coding"

    for exon in gene:
        exon_line.setAttribute("gene_biotype", gene_type)

    gene_line = GTF.Entry().fromGTF(gene[0])
    gene_line.feature = "gene"
    gene_line.start = gene_start
    gene_line.end = gene_end

    # Its not clear what the transcript_id for a gene line is.
    # Technically speaking, Ensembl GTFs do not follow the GTF
    # standard as GTF entries must include transcript_ids, but
    # the ensembl ones don't.

    if hasattr(gene_line, "gene_name"):
        gene_line.attributes["gene_name"] = gene[0].gene_name
    else:
        gene_line.attributes["gene_name"] = gene[0].gene_id

    gene_line.attributes["gene_biotype"] = gene_type

    print (str(gene_line)+"\n")
    for exon in gene:
        print(str(exon_line)+"\n")
</code></pre>

CC BY-SA 4.0
2018-12-11T10:33:29.077Z
That’s great to know. Since this pertains to using tools in R, it would help to have instructions to install the library and run for those unfamiliar with Python.
Tom Kelly ケリー・トム
CC BY-SA 4.0
2018-12-13T05:41:56.310Z
I've tested this script and made minor changes. This is what worked for me (to covert the GTF: I'm not sure if it's compatible with velocyto yet).
Tom Kelly ケリー・トム
CC BY-SA 4.0
2018-12-13T11:02:21.793Z
Thanks for the fixes. The things have changed recently and I lose track of what version things are on. I believe the CGAT team is working on simplifying the installation. A couple of points: I'm not quite sure what you are trying to achieve with the `if hasattr(gene_line, "gene_name)` bit: presumably if it already has a `gene_name` attribute, you don't need to set one? Also note that your resulting lines arn't going to have a `transcript_id`. While this is also true of Ensembl GTFs, it does mean these GTFs don't mean the  GTF specification.
Ian Sudbery
CC BY-SA 4.0
2018-12-13T11:12:00.660Z
I had issues with missing gene names and often TAIR IDs are used instead of gene symbols. So of this might be specific to Arabidopsis. Most important is the fixes to the loop for exons as these lines weren’t defined.
Tom Kelly ケリー・トム
CC BY-SA 4.0
2018-12-13T15:12:31.337Z
Its not that I don't understand that you might want to add a `gene_name` to the attributes, its just that currently you test for the presence of a `gene_name` and then set it to its already currently defined value if it already has one. This seems redundant.
Ian Sudbery
CC BY-SA 4.0
2018-12-11T10:23:37.900Z
2018-12-13T10:59:15.307Z
2018-12-13T10:59:15.307Z
Ian Sudbery
<h3>Why I need a compatible file</h3>

<p>I’m trying to run 
<a href="http://velocyto.org/" rel="nofollow noreferrer">velocyto</a> with the <a href="https://github.com/velocyto-team/velocyto.R" rel="nofollow noreferrer">R package</a> to analyse RNA velocity (cell trajectories) with single cell RNASeq data. I have performed single cell analysis from 10x Genomics data using cellranger.</p>

<p>I have successfully aligned the reads to get loom files and imported these into R. I can get the velocity from these files by following <a href="http://pklab.med.harvard.edu/velocyto/notebooks/R/chromaffin2.nb.html" rel="nofollow noreferrer">the vignettes</a>. However, I cannot reproduce the RNA velocity analysis based on “gene structure”.</p>

<p>I’m working with a different organism to the example (not human or mouse) so the annotation data provided in the example does not work. I have a GTF file for the latest annotation of this organism. However, it only contains “exon” and “CDS” as features. This appears to be the source of the problem. The “find.ip.sites” function requires a GTF with “features” = “gene” and one of the “attributes” to be “protein_coding”. These requirements are hard-coded into the velocyto.R function.</p>

<p>I have the following GTF files from the <a href="https://ics.hutton.ac.uk/atRTD/" rel="nofollow noreferrer">AtRTD2</a> dataset. The chromosome labels match the loom files in R.</p>

<pre><code>Chr1    TAIR10  exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    Araport11   exon    6788    7069    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7157    7450    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7564    7649    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7762    7835    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7942    7987    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8236    8325    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8417    8464    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8571    8737    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7315    7450    .   -   1   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7564    7649    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7762    7835    .   -   2   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7942    7987    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 8236    8325    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
</code></pre>

<p>The specifications for a <a href="https://asia.ensembl.org/info/website/upload/gff.html" rel="nofollow noreferrer">GTF or GFF2</a> are:</p>

<blockquote>
  <h3>Fields</h3>
  
  <p>Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'</p>
  
  <ol>
  <li><p>seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.</p></li>
  <li><p>source - name of the program that generated this feature, or the data source (database or project name)</p></li>
  <li><p>feature - feature type name, e.g. Gene, Variation, Similarity</p></li>
  <li><p>start - Start position of the feature, with sequence numbering starting at 1.</p></li>
  <li><p>end - End position of the feature, with sequence numbering starting at 1.</p></li>
  <li><p>score - A floating point value.</p></li>
  <li><p>strand - defined as + (forward) or - (reverse).</p></li>
  <li><p>frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..</p></li>
  <li><p>attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.</p></li>
  </ol>
</blockquote>

<h3>What I’m looking for</h3>

<p>Is there a way to derive a compatible GTF file containing protein coding genes from the exons and CDS? I want to produce a GTF which contains genes in the features and protein_coding in the attributes. Is it possible to do this with existing tools or scripts? </p>

<hr>

<h3>What I’ve tried so far</h3>

<p>I can modify the <a href="https://github.com/velocyto-team/velocyto.R/blob/666e1db1a6bebf042114eb9eb4c58f54cb002d45/R/momentum_routines.R#L2157-L2242" rel="nofollow noreferrer">source code</a> of the “find.ip.sites” function to run on my GTF file with these features missing. However, this requires running internal functions from the package written in Rcpp and means my workflow override future updates to the package. Running the rest of the vignette returns errors as no introns long enough have been identified (setting the thresholds lower is also incompatible with the general linear models called). Therefore I think it is better to generate a compatible GTF or GFF3 file rather than alerting the source code of the functions.</p>

<p>While intended for GTF files, the package functions can import GFF3 files, despite being described for GTF in the documentation. I’ve tried generating a GFF3 file using gffread from Cufflinks and replacing “feature” = “mRNA” with “gene” and adding “protein_coding”. This also returns errors when running the velocity algorithm. It does not work for either the GTF file used as an input for cellranger or the file generated by it. Is there a way to annotate protein coding genes and intron/exon boundaries based on a GTF file containing only exon and CDS annotations?</p>

<p>With cufflinks 2.2.1, a GFF3 was generated</p>

<pre><code>gffread - E AtRTD2_19April2016.gtf -o- &gt; AtRTD2_19April2016.gff
sed -i '/mRNA/s/gene_name=AT/gene_type=protein_coding;gene_name=AT/g' AtRTD2_19April2016.gff'
sed -i 's/mRNA/gene/g' AtRTD2_19April2016.gff
</code></pre>

<p>This is the GFF3 files that I've tried:</p>

<pre><code># gffread - E AtRTD2_19April2016.gtf -o- 
##gff-version3
Chr1    TAIR10  gene    3631    5899    .   +   .   ID=AT1G01010.1;geneID=AT1G01010;gene_type=protein_coding;gene_name=AT1G01010
Chr1    TAIR10  exon    3631    3913    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    3996    4276    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4486    4605    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4706    5095    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5174    5326    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5439    5899    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3760    3913    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3996    4276    .   +   2   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4486    4605    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4706    5095    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5174    5326    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5439    5630    .   +   0   Parent=AT1G01010.1
Chr1    Araport11   gene    6788    8737    .   +   .   ID=AT1G01020_P2;geneID=AT1G01020;gene_type=protein_coding;gene_name=AT1G01020    
Chr1    Araport11   exon    6788    7069    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7157    7450    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7564    7649    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7762    7835    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7942    7987    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8236    8325    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8417    8464    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8571    8737    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7315    7450    .   -   1   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7564    7649    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7762    7835    .   -   2   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7942    7987    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 8236    8325    .   -   0   Parent=AT1G01020_P2
</code></pre>

<p>The specifications for a <a href="https://asia.ensembl.org/info/website/upload/gff.html" rel="nofollow noreferrer">GTF or GFF2</a> are:     </p>

<p>The specifications for a <a href="https://asia.ensembl.org/info/website/upload/gff3.html" rel="nofollow noreferrer">GF3</a> are:</p>

<blockquote>
  <h3>Fields</h3>
  
  <p>Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'</p>
  
  <ol>
  <li><p>seqid - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seq ID must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.</p></li>
  <li><p>source - name of the program that generated this feature, or the data source (database or project name)</p></li>
  <li><p>type - type of feature. Must be a term or accession from the SOFA sequence ontology</p></li>
  <li><p>start - Start position of the feature, with sequence numbering starting at 1.</p></li>
  <li><p>end - End position of the feature, with sequence numbering starting at 1.</p></li>
  <li><p>score - A floating point value.</p></li>
  <li><p>strand - defined as + (forward) or - (reverse).</p></li>
  <li><p>phase - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..</p></li>
  <li><p>attributes - A semicolon-separated list of tag-value pairs, providing additional information about each feature. Some of these tags are predefined, e.g. ID, Name, Alias, Parent - see the GFF documentation for more details.</p></li>
  </ol>
</blockquote>

CC BY-SA 4.0
2018-12-10T15:15:25.017Z
It would probably help of you could [edit] your question and show us a few lines of each of the two GTF files so we can see what you mean more clearly.
terdon
CC BY-SA 4.0
2018-12-10T15:39:22.460Z
You probably can use AWK for this, is the "gene" feature similar as "CDS" or is it from "first exon start" to "last exon end".
benn
CC BY-SA 4.0
2018-12-10T18:24:23.847Z
I think the question and the "What I'm looking for" are a little confusing. Specifically, it sounds like you want to generate a GTF in "What I'm looking for" but you want to get the amino acid sequence of proteins from a GTF in the title. Is there any way you can clarify or make the title and "What I'm looking for" consistent?
conchoecia
CC BY-SA 4.0
2018-12-11T08:49:01.487Z
We have a [complete GTF for Arabidopsis in Ensembl Plants](ftp://ftp.ensemblgenomes.org/pub/plants/current/gtf/arabidopsis_thaliana/).
Emily_Ensembl
CC BY-SA 4.0
2018-12-18T01:51:46.953Z
Thanks @Emily_Ensembl, I was able to get velocyto to compute "global" RNA velocity based on this reference sequence. In case anyone else runs into this issue, the loom file needs to be re-computed with same GTF file with the shell utility to be compatible with the R commands.
Tom Kelly ケリー・トム
<scrnaseq><sequence-annotation><single-cell><gtf><10x-genomics>
Derive a GTF containing protein coding genes from a GTF file with Exons and CDS
# Derive a GTF containing protein coding genes from a GTF file with Exons and CDS

(Asked by: username_0 on 2018-12-10)

Why I need a compatible file
I’m trying to run 
velocyto with the R package to analyse RNA velocity (cell trajectories) with single cell RNASeq data. I have performed single cell analysis from 10x Genomics data using cellranger.
I have successfully aligned the reads to get loom files and imported these into R. I can get the velocity from these files by following the vignettes. However, I cannot reproduce the RNA velocity analysis based on “gene structure”.
I’m working with a different organism to the example (not human or mouse) so the annotation data provided in the example does not work. I have a GTF file for the latest annotation of this organism. However, it only contains “exon” and “CDS” as features. This appears to be the source of the problem. The “find.ip.sites” function requires a GTF with “features” = “gene” and one of the “attributes” to be “protein_coding”. These requirements are hard-coded into the velocyto.R function.
I have the following GTF files from the AtRTD2 dataset. The chromosome labels match the loom files in R.
Chr1    TAIR10  exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    Araport11   exon    6788    7069    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7157    7450    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7564    7649    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7762    7835    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7942    7987    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8236    8325    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8417    8464    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8571    8737    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7315    7450    .   -   1   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7564    7649    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7762    7835    .   -   2   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7942    7987    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 8236    8325    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";

The specifications for a GTF or GFF2 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.


What I’m looking for
Is there a way to derive a compatible GTF file containing protein coding genes from the exons and CDS? I want to produce a GTF which contains genes in the features and protein_coding in the attributes. Is it possible to do this with existing tools or scripts? 

What I’ve tried so far
I can modify the source code of the “find.ip.sites” function to run on my GTF file with these features missing. However, this requires running internal functions from the package written in Rcpp and means my workflow override future updates to the package. Running the rest of the vignette returns errors as no introns long enough have been identified (setting the thresholds lower is also incompatible with the general linear models called). Therefore I think it is better to generate a compatible GTF or GFF3 file rather than alerting the source code of the functions.
While intended for GTF files, the package functions can import GFF3 files, despite being described for GTF in the documentation. I’ve tried generating a GFF3 file using gffread from Cufflinks and replacing “feature” = “mRNA” with “gene” and adding “protein_coding”. This also returns errors when running the velocity algorithm. It does not work for either the GTF file used as an input for cellranger or the file generated by it. Is there a way to annotate protein coding genes and intron/exon boundaries based on a GTF file containing only exon and CDS annotations?
With cufflinks 2.2.1, a GFF3 was generated
gffread - E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gff
sed -i '/mRNA/s/gene_name=AT/gene_type=protein_coding;gene_name=AT/g' AtRTD2_19April2016.gff'
sed -i 's/mRNA/gene/g' AtRTD2_19April2016.gff

This is the GFF3 files that I've tried:
# gffread - E AtRTD2_19April2016.gtf -o- 
##gff-version3
Chr1    TAIR10  gene    3631    5899    .   +   .   ID=AT1G01010.1;geneID=AT1G01010;gene_type=protein_coding;gene_name=AT1G01010
Chr1    TAIR10  exon    3631    3913    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    3996    4276    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4486    4605    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4706    5095    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5174    5326    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5439    5899    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3760    3913    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3996    4276    .   +   2   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4486    4605    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4706    5095    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5174    5326    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5439    5630    .   +   0   Parent=AT1G01010.1
Chr1    Araport11   gene    6788    8737    .   +   .   ID=AT1G01020_P2;geneID=AT1G01020;gene_type=protein_coding;gene_name=AT1G01020    
Chr1    Araport11   exon    6788    7069    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7157    7450    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7564    7649    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7762    7835    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7942    7987    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8236    8325    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8417    8464    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8571    8737    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7315    7450    .   -   1   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7564    7649    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7762    7835    .   -   2   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7942    7987    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 8236    8325    .   -   0   Parent=AT1G01020_P2

The specifications for a GTF or GFF2 are:     
The specifications for a GF3 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqid - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seq ID must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
type - type of feature. Must be a term or accession from the SOFA sequence ontology
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
phase - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attributes - A semicolon-separated list of tag-value pairs, providing additional information about each feature. Some of these tags are predefined, e.g. ID, Name, Alias, Parent - see the GFF documentation for more details.

**Tags:** scrnaseq, sequence-annotation, single-cell, gtf, 10x-genomics
> **Comments:**
> **username_1**: It would probably help of you could [edit] your question and show us a few lines of each of the two GTF files so we can see what you mean more clearly.  
> **username_2**: You probably can use AWK for this, is the "gene" feature similar as "CDS" or is it from "first exon start" to "last exon end".  
> **username_3**: I think the question and the "What I'm looking for" are a little confusing. Specifically, it sounds like you want to generate a GTF in "What I'm looking for" but you want to get the amino acid sequence of proteins from a GTF in the title. Is there any way you can clarify or make the title and "What I'm looking for" consistent?  
> **username_4**: We have a [complete GTF for Arabidopsis in Ensembl Plants](ftp://ftp.ensemblgenomes.org/pub/plants/current/gtf/arabidopsis_thaliana/).  
> **username_0**: Thanks @Emily_Ensembl, I was able to get velocyto to compute "global" RNA velocity based on this reference sequence. In case anyone else runs into this issue, the loom file needs to be re-computed with same GTF file with the shell utility to be compatible with the R commands.  

---

## 1 answers

### Answer by username_5 on 2018-12-11. Score: 4 | Accepted answer

In your case, I would definitely suggest following @Emily_Ensembl's advice and using the Arabidopsis GTF from Ensembl. But for future reference, if an Ensembl GTF wasn't available, you could build something like this using the gtf class from cgat
The cgat module and dependancies can be installed by following the instructions here. Specifically, you can download and run their installation script. This will set up a dedicated conda environment to run the versions of dependancies needed. Further documentation can be found here.
# download installation script:
curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/master/install.sh

# install the development version (recommended, no production version yet):
bash install.sh --devel

This requires python3 and anaconda to install. You may need to run this script in the conda environment:
# show available environments
conda info --envs
# activate by path
source activate /home/user/local/bin/conda-install/envs/cgat-a

Then run the follow as a python3 script. Copy the contents in a file such as convert_gtf.py. Then run it in the terminal on your gtf file:
python convert_gtf.py genes.gtf > genes_new.gtf

#python 3.6.4

import sys

from cgat import GTF
from cgatcore import iotools

for gene in GTF.flat_gene_iterator(
     GTF.iterator(iotools.open_file(sys.argv[1])), strict = False):
    gene_start = min (e.start for e in gene)
    gene_end = max(e.end for e in gene)

    if any(e.feature == "CDS" for e in gene):
        gene_type = "protein_coding"
    else:
        gene_type = "non_coding"

    for exon in gene:
        exon_line.setAttribute("gene_biotype", gene_type)

    gene_line = GTF.Entry().fromGTF(gene[0])
    gene_line.feature = "gene"
    gene_line.start = gene_start
    gene_line.end = gene_end

    # Its not clear what the transcript_id for a gene line is.
    # Technically speaking, Ensembl GTFs do not follow the GTF
    # standard as GTF entries must include transcript_ids, but
    # the ensembl ones don't.

    if hasattr(gene_line, "gene_name"):
        gene_line.attributes["gene_name"] = gene[0].gene_name
    else:
        gene_line.attributes["gene_name"] = gene[0].gene_id

    gene_line.attributes["gene_biotype"] = gene_type

    print (str(gene_line)+"\n")
    for exon in gene:
        print(str(exon_line)+"\n")

> **Comments:**
> **username_0**: That’s great to know. Since this pertains to using tools in R, it would help to have instructions to install the library and run for those unfamiliar with Python.  
> **username_0**: I've tested this script and made minor changes. This is what worked for me (to covert the GTF: I'm not sure if it's compatible with velocyto yet).  
> **username_5**: Thanks for the fixes. The things have changed recently and I lose track of what version things are on. I believe the CGAT team is working on simplifying the installation. A couple of points: I'm not quite sure what you are trying to achieve with the `if hasattr(gene_line, "gene_name)` bit: presumably if it already has a `gene_name` attribute, you don't need to set one? Also note that your resulting lines arn't going to have a `transcript_id`. While this is also true of Ensembl GTFs, it does mean these GTFs don't mean the  GTF specification.  
> **username_0**: I had issues with missing gene names and often TAIR IDs are used instead of gene symbols. So of this might be specific to Arabidopsis. Most important is the fixes to the loop for exons as these lines weren’t defined.  
> **username_5**: Its not that I don't understand that you might want to add a `gene_name` to the attributes, its just that currently you test for the presence of a `gene_name` and then set it to its already currently defined value if it already has one. This seems redundant.  

---

