Hi,

I would like to use edgeR to estimate the RPKM values.  I have (1) read counts files estimated by HTSeq-count, and (2) a transcript length file. I know how to estimate CPM in edgeR, using below command lines. However, I don’t know how to estimate RPKM values based on the files I have. Could you show me how to run the command lines in edgeR? Many thanks.

Best,

Gary

Transcript count files

http://68.181.92.180/~Gary/temporal/T8_Early_R_T.txt

http://68.181.92.180/~Gary/temporal/T8_Late_R_T.txt

http://68.181.92.180/~Gary/temporal/T9_Early_R_T.txt

http://68.181.92.180/~Gary/temporal/T9_Late_R_T.txt

Transcript length file

http://68.181.92.180/~Gary/temporal/hg19refGene20160206_TranscriptLength.txt

> targets=readTargets()

> targets

             files group  description

1 T8_Early_R_T.txt Early T8_Early_R_T

2  T8_Late_R_T.txt  Late  T8_Late_R_T

3 T9_Early_R_T.txt Early T9_Early_R_T

4  T9_Late_R_T.txt  Late  T9_Late_R_T

> d=readDGE(targets,skip=0,comment.char="!")

> keep=rowSums(cpm(d)>1)>=2

> d=d[keep,,keep.lib.sizes=FALSE]

> d=calcNormFactors(d,method="TMM")

> dim(d)

[1] 8202    4

> detags=rownames(d)

> detagsfile=cpm(d)[detags,]

> write.table(detagsfile,file="CPM_MSC_8202ene.txt")

> rpkm<-rpkm(d)

Error in rpkm.DGEList(d) : Gene lengths not found

> rpkm<-rpkm(d,gene.lengths)

Error in rpkm.DGEList(d, gene.lengths) : object 'gene.lengths' not found

Well, it would be a good start to load the gene lengths into memory.

gene.data <- read.table("hg19refGene20160206_TranscriptLength.txt", header=TRUE)

You then need to match the gene names in d to the ordering of gene.data:

m <- match(rownames(d), gene.data$Transcript)
gene.lengths <- gene.data$TranscriptLength[m]

... and then you should be able to run rpkm on that.