SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
[DEXSeq] problem with estimateSizeFactors yerbol Bioinformatics 15 11-08-2013 01:57 PM
Problem with DEXSeq read.HTSeqCounts marianaboroni Bioinformatics 0 03-13-2013 11:43 AM
question with read.HTSeqCounts function in DEXSeq package bioliyezhang Bioinformatics 2 12-11-2012 06:33 AM
Problem working with DEXSeq package senkewiczs Bioinformatics 0 07-17-2012 08:01 AM
using DEXSeq one gene with only one exon - problem Jane M RNA Sequencing 1 08-04-2011 06:09 AM

Reply
 
Thread Tools
Old 01-09-2014, 02:43 AM   #1
Erwan
Junior Member
 
Location: Nantes, France

Join Date: Jan 2014
Posts: 8
Default DEXSeq read.HTSeqCounts problem

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 :

Code:
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:

Code:
_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
Erwan is offline   Reply With Quote
Old 01-10-2014, 12:25 AM   #2
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

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
areyes is offline   Reply With Quote
Old 01-10-2014, 06:33 AM   #3
Erwan
Junior Member
 
Location: Nantes, France

Join Date: Jan 2014
Posts: 8
Default

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.

Code:
 ## 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()
Code:
> 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 :


Code:
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

Code:
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,
Attached Files
File Type: gz run_DEXseq.R.tar.gz (1,017 Bytes, 5 views)
Erwan is offline   Reply With Quote
Old 01-13-2014, 01:16 AM   #4
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

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
areyes is offline   Reply With Quote
Old 01-15-2014, 12:24 AM   #5
Erwan
Junior Member
 
Location: Nantes, France

Join Date: Jan 2014
Posts: 8
Default

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,
Erwan is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:37 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO