5
$\begingroup$

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 '.'

  1. 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.

  2. source - name of the program that generated this feature, or the data source (database or project name)

  3. feature - feature type name, e.g. Gene, Variation, Similarity

  4. start - Start position of the feature, with sequence numbering starting at 1.

  5. end - End position of the feature, with sequence numbering starting at 1.

  6. score - A floating point value.

  7. strand - defined as + (forward) or - (reverse).

  8. 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..

  9. 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 '.'

  1. 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.

  2. source - name of the program that generated this feature, or the data source (database or project name)

  3. type - type of feature. Must be a term or accession from the SOFA sequence ontology

  4. start - Start position of the feature, with sequence numbering starting at 1.

  5. end - End position of the feature, with sequence numbering starting at 1.

  6. score - A floating point value.

  7. strand - defined as + (forward) or - (reverse).

  8. 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..

  9. 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.

Improve this question
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – terdon on strike ♦
    Dec 10, 2018 at 15:15
  • $\begingroup$ 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". $\endgroup$
    – benn
    Dec 10, 2018 at 15:39
  • $\begingroup$ 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? $\endgroup$
    – conchoecia
    Dec 10, 2018 at 18:24
  • 3
    $\begingroup$ We have a complete GTF for Arabidopsis in Ensembl Plants. $\endgroup$
    – Emily_Ensembl
    Dec 11, 2018 at 8:49
  • $\begingroup$ 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. $\endgroup$
    – Tom Kelly
    Dec 18, 2018 at 1:51

1 Answer 1

Reset to default
4
$\begingroup$

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")
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Tom Kelly
    Dec 11, 2018 at 10:33
  • $\begingroup$ 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). $\endgroup$
    – Tom Kelly
    Dec 13, 2018 at 5:41
  • $\begingroup$ 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. $\endgroup$
    – Ian Sudbery
    Dec 13, 2018 at 11:02
  • $\begingroup$ 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. $\endgroup$
    – Tom Kelly
    Dec 13, 2018 at 11:12
  • $\begingroup$ 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. $\endgroup$
    – Ian Sudbery
    Dec 13, 2018 at 15:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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