get gene lines from gtf file

Question

I would like to retrieve gene lines from a GTF file for which I only have exons & transcripts lines (output from Cufflinks) and alternative splicing possible. I need gene lines for compatibility with a pipeline dealing with GTF in Ensembl format.

Ideally, I would have gene lines representing the longest transcript - but I am open to discussion about that. The aim would be to separate the genome in genic & intergenic portions - boundaries being defined by start-end coordinates of the genes.

Sample input (e.g. gene "CUFF.105" with 5 transcripts):

cat subset_105.test | awk '$3=="transcript"{print $0}'

tig00000046 Cufflinks   transcript  26170   42766   202 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.3"; FPKM "0.2320304094"; frac "0.197788"; conf_lo "0.164885"; conf_hi "0.299176"; cov "5.188855";
tig00000046 Cufflinks   transcript  26170   42766   266 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.1"; FPKM "0.3061260755"; frac "0.249779"; conf_lo "0.228860"; conf_hi "0.383392"; cov "6.845843";
tig00000046 Cufflinks   transcript  26170   39469   470 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.2"; FPKM "0.5403628578"; frac "0.161731"; conf_lo "0.372512"; conf_hi "0.708214"; cov "12.084038";
tig00000046 Cufflinks   transcript  28928   39469   1000    +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.4"; FPKM "1.1485378983"; frac "0.233911"; conf_lo "0.853387"; conf_hi "1.443688"; cov "25.684547";
tig00000046 Cufflinks   transcript  29614   42766   181 +   .   gene_id "CUFF.105"; transcript_id "CUFF.105.5"; FPKM "0.2087076495"; frac "0.156792"; conf_lo "0.137112"; conf_hi "0.280304"; cov "4.667292";

What it came to my mind is:

• sort the file by $4 (numeric: n) and $5 (reverse: nr) (sort -k4,4n -k5,5nr)
• take the 1st line
• remove unnecessary fields from $9

But this would fail for transcripts annotated in the reverse strand, and I won't be able to deal (how to fill) the $6 (score - a floating point value).

Any ideas of an existing method that could help me retrieve those gene lines?

Options to retrieve those lines directly from Cufflinks would also be accepted (I didn't find a suitable parameter in the manual).

Answer 1

Assuming none of the genes have duplicate entries across chromosomes (or strand or other biologically implausible things) in R one can:

library(rtracklayer)
library(GenomicRanges)
gtf = import.gff("file.gtf")
grl = split(gtf, gtf$gene_id)  # This produces a GRangesList
grl = endoapply(grl, function(x) {
    foo = x[1]
    foo$type = 'gene'
    start(foo) = min(start(x))
    end(foo) = max(end(x))
    return(c(foo, x))
    })
gr = unlist(gr)

You can then either write gr out to a file (either with export() or just printing the rows).