DEXSEQ Prepare Annotation File and R output

[gokhulkrishnakilaru]

Hi Folks,

I have downloaded the DEXSEQ package from Bioconductor.

When I try to make the annotation file using the dexseq_prepare_annotations.py script, a gff file is generated but is of zero KB in size.

I tried with the following files

Mus_musculus.NCBIM37.64.gtf
Mouse_UCSC_Refgene.gtf (Refgene from UCSC)
Mouse_UCSC_RefFlat.gtf (Refflat from UCSC)

And, all I got was a zero KB file.

The command I tried was

python dexseq_prepare_annotation.py Mus_musculus.NCBIM37.64.gtf Mus_musculus.NCBIM37.64.gff

The above command works for Mus_musculus.GRCm38.68.gtf downloaded from Ensemble and I don't understand why it wouldn't work for the NCBI gtf from the same Ensemble.

Also, the R output from running DEXSEQ using Mus_musculus.GRCm38.68.gf and one wild type and one knock out sample (no replicates) didn't give me anything. All I see is NA in all the columns. I am using the bam files from Tophat aligned to Refseq annotation.

What am I doing wrong? Any pointers are highly appreciated.

[areyes]

Have you checked the size of your first file? Looks like you are replacing your input file with the output file.

Could you please include a reproducible code for your R code? with the output of the sessionInput()? Also, I would not have many hopes in the results without replicates

Alejandro Reyes

[gokhulkrishnakilaru]

Have you checked the size of your first file? Looks like you are replacing your input file with the output file.

Could you please include a reproducible code for your R code? with the output of the sessionInput()? Also, I would not have many hopes in the results without replicates

Alejandro Reyes

Hi Alejandro,

Yes, I did check the size of the input file. I am changing the extension. The input file has GTF and the output has GFF as its extension.

My R code is as follows

library(DEXSeq)
options(digits=3)
setwd("/test/dexseq/")
library(DEXSeq)
rm(list=ls())
annotationfile = file.path("/test/dexseq/Mus_musculus.GRCm38.68.gff")
annotationfile
samples = data.frame(condition = c("WT", "KO"),replicate=c(1,1),row.names=c("WildType", "KnockOut"),stringsAsFactors=TRUE,check.names = FALSE)
samples
fullFilenames<- list.files("/test/dexseq/",full.names=TRUE,pattern="DEXSEQ.txt")
fullFilenames
ecs<- read.HTSeqCounts(countfiles = fullFilenames,design = samples,flattenedfile = annotationfile)
head(counts(ecs))
head(fData(ecs))

All I see is NA and the estimate size factor is also giving out NA.

[areyes]

ups, my bad in the gtf extensions thing...

do your files contain NAs also?

[gokhulkrishnakilaru]

ups, my bad in the gtf extensions thing...

does your files contain NAs also?

No. My files have either a value or 0 for nothing. Also, I see another error saying

error in scan(file what nmax sep dec quote skip nlines na.strings line 1 did not have 3 elements

I looked the tail of my counts file and it has got some four lines in the last saying _ambiguous, _lowqual etc.

I deleted those lines and it gives me another error saying

Error in round(countData) : Non-numeric argument to mathematical function

.

Any pointers to these issues. This is my counts file's head

"ENSMUSG00000000001"    :001	1
"ENSMUSG00000000001"    :002	0
"ENSMUSG00000000001"    :003	0
"ENSMUSG00000000001"    :004	1
"ENSMUSG00000000001"    :005	0
"ENSMUSG00000000001"    :006	0
"ENSMUSG00000000001"    :007	0
"ENSMUSG00000000001"    :008	0

[areyes]

I see, I think the files you are using as input are causing some problems with the output of our htseq python scripts. I will check what is going on. In the meantime you can reformat your files to look more like this:

FBgn0000003:001	0
FBgn0000008:001	0
FBgn0000008:002	0
FBgn0000008:003	0
FBgn0000008:004	1
FBgn0000008:005	4
FBgn0000008:006	1
FBgn0000008:007	18
FBgn0000008:008	4
FBgn0000008:009	16

Then it should be fine!

[areyes]

By the way, where can I download the annotation files you used?

[gokhulkrishnakilaru]

I see, I think the files you are using as input are causing some problems with the output of our htseq python scripts. I will check what is going on. In the meantime you can reformat your files to look more like this:

FBgn0000003:001	0
FBgn0000008:001	0
FBgn0000008:002	0
FBgn0000008:003	0
FBgn0000008:004	1
FBgn0000008:005	4
FBgn0000008:006	1
FBgn0000008:007	18
FBgn0000008:008	4
FBgn0000008:009	16

Then it should be fine!

Perfect. That helps me. Also, what about those last four lines with the underscore sign and a numerical value. Can I delete them?

[gokhulkrishnakilaru]

By the way, where can I download the annotation files you used?

ftp://ftp.ensembl.org/pub/release-68/gtf/mus_musculus

That is where I got the one that worked for me.

You can use genome.ucsc.edu and go to tables section. Choose mouse and refseq genes and then refFlat or refGene. Select format to be GTF and if you are successful in preparing the annotations file. Please upload it somewhere or I can invite you to my dropbox. So, that way I have a refseq annotation file.

Thanks for the support, my friend.

[areyes]

You could, but they are also deleted automatically in the function "read.HTSeqCounts"!

[gokhulkrishnakilaru]

You could, but they are also deleted automatically in the function "read.HTSeqCounts"!

Hi Alejandro,

I was successful in making the counts file as you suggested. I ran the script. The following are my errors. Any pointers that could be of help?

ecs<- estimateSizeFactors(ecs)
> ecs<- estimateDispersions(ecs)
Dispersion estimation. (Progress report: one dot per 100 genes)
Error in FUN(c("ENSMUSG00000000078", "ENSMUSG00000000134", "ENSMUSG00000000182",  : 
  Underdetermined model; cannot estimate dispersions. Maybe replicates have not been properly specified.
In addition: Warning messages:
1: In .local(object, ...) :
  Exons with less than 11 counts will be discarded. For more details read the documentation, parameter minCount
2: In .local(object, ...) :
  Genes with more than 70 testable exons will be kicked out of the analysis. For more details read the documentation, parameter maxExon

I was looking at this link - http://seqanswers.com/forums/archive...p/t-21212.html. Can I delete that line for my case?

[gokhulkrishnakilaru]

Any thoughts anybody?

Sorry mods, for bumping up posts.

Urgent task. So, had to.

[areyes]

Hi gokhulkrishnakilaru,

The error talks by its own: "Underdetermined model; cannot estimate dispersions. Maybe replicates have not been properly specified.", you do not have replicates. Sorry that I can not help.

Alejandro

[gokhulkrishnakilaru]

Hi gokhulkrishnakilaru,

The error talks by its own: "Underdetermined model; cannot estimate dispersions. Maybe replicates have not been properly specified.", you do not have replicates. Sorry that I can not help.

Alejandro

Thanks Alejandro,

So no dexseq could work without replicates?

Is that the conclusion?

Is there a possibility to change the declaration while specifying the replicates in this section

samples = data.frame(condition = c("WT", "KO"),replicate=c(1,1),row.names=c("WildType", "KnockOut"),stringsAsFactors=TRUE,check.names = FALSE)

[areyes]

The motivation of the development of DESeq and DEXSeq is being able to estimate biological variability between replicates, and take this into account to call differentially expressed genes or exons. If you don´t have replicates, you do not know if the changes that you are observing are due to biological variation or due to the differences in your genotypes. In any experiment is crucial to do replicates, this is the only way to guarantee reproducibility on your differential expressed calls. For more details, you could check:

http://www.nature.com/nbt/journal/v29/n7/pdf/nbt.1910.pdf

the discussion of our DEXSeq paper:

http://genome.cshlp.org/content/22/10/2008.full.pdf+html