DEXSeq read.HTSeqCounts problem

[Erwan]

Hi,

First I'd like to say that I'm very new to bioinformatics so don't hesitate to ask me for further details if I'm not clear enough.

I'm trying to use DEXSeq for an analysis but I'm running into the following error :

Error in `rownames<-`(`*tmp*`, value = character(0)) : 
  attempt to set 'rownames' on an object with no dimensions
Calls: read.HTSeqCounts -> rownames<-
Execution halted

As inputs, I use a gene annotation file from UCSC and a BAM file produced by Tophat2.

The error occurs in the third phase of the DEXSeq-hts.sh script ("Differential testing"), during the execution of the R script "run_DEXseq.R" and especially the call of the read.HTSeqCounts function. The script has previously created a specific annotation file and a "read counting" file that is used in the R script.

Here are the first lines of the "read counting" file if it can be of any help:

_ambiguous      0
_empty  1922672
_lowaqual       0
_notaligned     0
chr10_1190007I07Rik-:001        0
chr10_1190007I07Rik-:002        0
chr10_1190007I07Rik-:003        0
chr10_1190007I07Rik-:004        0
chr10_1190007I07Rik-:005        0
chr10_1500009L16Rik+:001        0
chr10_1500009L16Rik+:002        0
chr10_1500009L16Rik+:003        1
chr10_1500009L16Rik+:004        1
chr10_1500009L16Rik+:005        0

If someone has already encountered this error, or has any clue about it, I'm more than interested !

Erwan

[areyes]

Hi Erwan,

Could you include a reproducible code that lead to the error? also could you include the output of your "sessionInfo()"?

What scripts are you referring when you mention "DEXSeq-hts.sh"?

Alejandro

[Erwan]

Hi Alejandro,

Yes I forgot to mention that I'm using Galaxy and the script DEXSeq-hts.sh is part of the DESeq package we can get from the Galaxy main tool shed, it automatizes the steps of a DEXSeq analysis.

But I've also tried to use DEXSeq in the standard way (as described in your document "Inferring differential exon usage in RNA-Seq data with the DEXSeq package") and I'm getting the same issue.

I also tried to use ENSEMBL files rather than UCSC files as it is recommended ==> same issue.

Another thing that may be useful. When I run the dexseq_count.py script I get a whole bunch of warnings like this : "/usr/local/lib/python2.7/dist-packages/HTSeq/__init__.py:598: UserWarning: Read HWI-1KL149:612C11ACXX:8:2316:21358:21144 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
"which could not be found. (Is the SAM file properly sorted?)" )"

And I even get these warnings when I sort the bam file with samtools.

Eventually, here is what you asked for. I only put a few lines of the R script run by DEXSes-hts.sh but if you want to see more, I've attached it to this post.

 ## header name
    test_pair_name <- c(paste(pw_tests[[i]][1], "__vs__", pw_tests[[i]][2], sep=""))
    print(test_pair_name)

    sub.data <- subset(condsTable, (conditions %in% c(pw_tests[[i]][1],pw_tests[[i]][2])))
    sub.data[[1]]<-as.factor(sub.data[[1]])

    ecs = read.HTSeqCounts(countfiles=file.path(EXTRAPATH, row.names(sub.data)), design=sub.data, flattenedfile=annodb)

The error comes at the execution of the read.HtSeqCounts function. Another thing that may be useful, the print(test_pair_name) instruction prints the following line : [1] "1__vs__NA"

Here is my sessionInfo()

> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
 [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] DEXSeq_1.8.0       Biobase_2.22.0     BiocGenerics_0.8.0

loaded via a namespace (and not attached):
 [1] biomaRt_2.18.0       Biostrings_2.30.1    bitops_1.0-6        
 [4] GenomicRanges_1.14.4 hwriter_1.3          IRanges_1.20.6      
 [7] RCurl_1.95-4.1       Rsamtools_1.14.2     statmod_1.4.18      
[10] stats4_3.0.2         stringr_0.6.2        XML_3.98-1.1        
[13] XVector_0.2.0        zlibbioc_1.8.0

Here are the last lines of my sam file :

HWI-1KL149:61:D2C11ACXX:8:2302:21370:41437      417     Y       2575578 3       101M    3       5860463 0       AGCCACTCAGTACTGAGAGTGGAGATAGAGGTCAAATCCCCAGTGTGCAAAAATAGGGGCATTACACAAAACCTTCTCCCAGGCTCGTCACCCAGAGGGGT     D@CCFFFFFGHFHHJIJCGIIGDE@GHIIJI<?EEHGHGIGGIHGGHBGHIFHIGHEGG68CHIGEIIJ>CEHFFFFFBECE@BBDDDDDBC<2<BB@BB#   AS:i:-10        XN:i:0    XM:i:2  XO:i:0  XG:i:0  NM:i:2  MD:Z:0G5C94     YT:Z:UU NH:i:2  HI:i:1
HWI-1KL149:61:D2C11ACXX:8:1303:15000:65350      129     Y       2622474 50      101M    4       147360561       0       ACTTCACCTTCAGCCAGCACCTCAGCAGGTCTTGGCTCTTTTCCTGGAGGATTGACCCATACCTTCATTCCTGATGGGTCTGCGTCCTTTGTCATCCTGCT     FCCCFFFFFHHHHHJJJJJJJJJJJJJJJHIJJJJJJIJJJJJJJJJIIJIIIJJJJJJJJJJJJJJJJJHHHHHHFFFFFEEDDDDDDDDDDEEDDDDDD   AS:i:0  XN:i:0    XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:101        YT:Z:UU NH:i:1

And here the last lines from the gff files produced by the prepare annotation python script

Y	dexseq_prepare_annotation.py	exonic_part	90838869	90839177	.	-	.	transcripts "ENSMUST00000179623"; exonic_part_number "001"; gene_id "ENSMUSG00000096850"
Y	dexseq_prepare_annotation.py	aggregate_gene	90843934	90844074	.	+	.	gene_id "ENSMUSG00000095326"
Y	dexseq_prepare_annotation.py	exonic_part	90843934	90844074	.	+	.	transcripts "ENSMUST00000178550"; exonic_part_number "001"; gene_id "ENSMUSG00000095326"

Erwan,

[areyes]

Hi Erwan,

I think the error is not coming from DEXSeq but from that R script used by Galaxy. Maybe the people that wrote that script are able to tell you more. My guess is that one of the conditions that you are specifying is NA, so the script is looking for non-existing files, I will change the DEXSeq code so that is more explicit if the user specifies files that does not exist.

Alejandro

[Erwan]

Hi Alejandro,

I eventually figured it out. The problem was not due to the Galaxy script but to me. It happens that I was only providing one input file to DEXSeq, that can be a little bit problematic for a differential tool... I still don't know how I could have missed that, next time I'll read the documentation more carefully.
Anyway, thanks for your help.

Erwan,