[dpryan]

"Error in letterFrequency(getSeq(FASTA, reducedGTF), "GC") :
error in evaluating the argument 'x' in selecting a method for function 'letterFrequency': Error in value[[3L]](cond) :
record 1666 (chr6_ssto_hap7:1871448-1871615) failed
file: genome.fa
Calls: getSeq ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous>
Execution halted
"

GTF files match the "top assembly", rather than the "primary assembly". Either download the matching fasta file or remove the haplotype patches from your GTF.

[super0925]

GTF files match the "top assembly", rather than the "primary assembly". Either download the matching fasta file or remove the haplotype patches from your GTF.

Sorry D, I didn't totally get what you said.
The genes.gtf and genome.fa are downloaded from the latest UCSC homo sapiens hg19 genome.
How could I 'download the matching fasta' or 'remove the haplotype patches' to use your script?
Thank you!

PS: I only change your script to
GTF <- import.gff(GTFfile, format="gtf", genome="hg19", asRangedData=F, feature.type="exon")
grl <- reduce(split(GTF, elementMetadata(GTF)$gene_id))
reducedGTF <- unlist(grl, use.names=T)

[dpryan]

If you just need the lengths then remove the code that calculates GC content. Something like the following should work:

#!/usr/bin/Rscript
library(GenomicRanges)
library(rtracklayer)
library(Rsamtools)

GTFfile = "something.GTF"

#Load the annotation and reduce it
GTF <- import.gff(GTFfile, format="gtf", genome="GRCm38.71", asRangedData=F, feature.type="exon")
grl <- reduce(split(GTF, elementMetadata(GTF)$gene_id))
reducedGTF <- unlist(grl, use.names=T)
elementMetadata(reducedGTF)$gene_id <- rep(names(grl), elementLengths(grl))

elementMetadata(reducedGTF)$widths <- width(reducedGTF)

#Create a list of the ensembl_id/GC
calc_GC_length <- function(x) {
    sum(elementMetadata(x)$widths)
}
output <- sapply(split(reducedGTF, elementMetadata(reducedGTF)$gene_id), calc_GC_length)

[super0925]

If you just need the lengths then remove the code that calculates GC content. Something like the following should work:

#!/usr/bin/Rscript
library(GenomicRanges)
library(rtracklayer)
library(Rsamtools)

GTFfile = "something.GTF"

#Load the annotation and reduce it
GTF <- import.gff(GTFfile, format="gtf", genome="GRCm38.71", asRangedData=F, feature.type="exon")
grl <- reduce(split(GTF, elementMetadata(GTF)$gene_id))
reducedGTF <- unlist(grl, use.names=T)
elementMetadata(reducedGTF)$gene_id <- rep(names(grl), elementLengths(grl))

elementMetadata(reducedGTF)$widths <- width(reducedGTF)

#Create a list of the ensembl_id/GC
calc_GC_length <- function(x) {
    sum(elementMetadata(x)$widths)
}
output <- sapply(split(reducedGTF, elementMetadata(reducedGTF)$gene_id), calc_GC_length)

Hi D
Brilliant! It works.
I have just changed one command
GTF <- import.gff(GTFfile, format="gtf", genome="hg19", asRangedData=F, feature.type="exon")
However, the result only contain 25369 genes.
However, there are 23170 rows in the count table generated by HTseq, which means there are 23170 genes in the human genome.
Why does the dimension are different? I don't think rpkm() could go further if the dimension are different...
Thank you!

[dpryan]

There are some entries that htseq-count will skip (those without exons?), which presumably is why this happens. Anyway, you need to ensure that everything is in the same order anyway, so make sure to write meaningful row names and use those with match().

[gringer]

Convert your BAM file to FASTQ, then use Kallisto.

The crazy thing is that this might actually be quicker than other methods.

[super0925]

There are some entries that htseq-count will skip (those without exons?), which presumably is why this happens. Anyway, you need to ensure that everything is in the same order anyway, so make sure to write meaningful row names and use those with match().

Thank you D. Very helpful.
And I also found Cufflinks could calculate RPKM directly.
by
cufflinks -o OutDir -g hg19/genes.gtf TopHat/accepted_hits.bam

[LacquerHead]

Hi

So the script worked for me and I got an mm9.fa.fai file as output, wondering how this can be used with edgeR to generate the RPKMs? Thanks.