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.