DEXseq errors

[Anomilie]

I've come across two errors while running DEXSeq that I haven't been able to resolve.

Firstly, when looking at the results using DEXSeqResults, I find some NA in the expression value and log2fold columns, while p.value and padj are available.

pvalue padj MT WT log2fold_MT_WT
<numeric> <numeric> <numeric> <numeric> <numeric>
2010010A06Rik:1 1.274581e-07 2.968773e-05 NA NA NA
2010010A06Rik:2 2.271527e-04 1.387325e-02 NA NA NA

Secondly, I have not been able to generate the HTML reports using the DEXSeqHTML command. It produces the following error:

DEXSeqHTML(dxr1, fitExpToVar="condition", FDR=0.05,color=c("#FF000080", "#0000FF80"),path="DEX_seq/",file="DEU_res_SO.html")
Error in data.frame(..., check.names = FALSE) :
arguments imply differing number of rows: 214419, 0

Has anyone come across something similar or have any suggestions on how to resolve these issues?

## makeTranscriptDbFromGFF
gffFile <- makeTranscriptDbFromGFF("Genome_files/Mus_musculus/UCSC/mm10/Annotation/Genes/genes.gtf", format="gtf")

## preparing exonic parts
exonicParts <- disjointExons(gffFile, by="exon", aggregateGenes=TRUE)


align <- "Samples"

files <- list.files(path=align, pattern="*.bam", full.names=T, recursive=FALSE)

bf1 <- BamFileList(c(files),index=character())


genehits <- summarizeOverlaps(exonicParts, bf1, mode="Union", ignore.strand=FALSE, singleEnd=FALSE, inter.feature=FALSE, fragments=TRUE)

colData(genehits)$condition <- c("WT", "MT", "WT", "MT", "WT", "MT")
raw_dat <- assays(genehits)$counts

conds <- c("WT1", "MT1", "WT2", "MT2", "WT3", "MT3")
colnames(raw_dat) <- conds
## reorder columns
raw_data <- cbind(raw_dat[,c(1,3,5)], raw_dat[,c(2,4,6)]) 


geneID <- exonicParts$gene_id #-> contains gene id
g <- unlist(geneID)
exonID <- exonicParts$exonic_part # contains exon number

nam <- paste(g,exonID, sep=":")

rownames(raw_dat) <- nam

#############################
### Generate DEXSeq object ##
#############################

sampleTable <- data.frame(row.names = conds, condition= factor(c("WT", "WT", "WT", "MT", "MT", "MT")))
design <- formula(~ sample + exon + condition:exon)
dxd <- DEXSeqDataSet(raw_data,sampleTable, design, featureID= as.character(exonID), groupID= g)

######################
### Perform testing ##
######################

dxd <- estimateSizeFactors(dxd)
dxd <- estimateDispersions(dxd)

dxd <- testForDEU(dxd)
dxd <- estimateExonFoldChanges(dxd, fitExpToVar="condition")
dxr1 = DEXSeqResults( dxd )


##################
## save results ##
##################

DEXSeqHTML(dxr1, fitExpToVar="condition", FDR=0.05,color=c("#FF000080", "#0000FF80"),path="/DEX_seq/",file="DEU_res_SO.html")

My sessionInfo is as following

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-unknown-linux-gnu (64-bit)

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

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

other attached packages:
 [1] GenomicAlignments_1.0.1   BSgenome_1.32.0          
 [3] Rsamtools_1.16.0          Biostrings_2.32.0        
 [5] XVector_0.4.0             biomaRt_2.20.0           
 [7] rtracklayer_1.24.2        GenomicFeatures_1.16.2   
 [9] AnnotationDbi_1.26.0      DEXSeq_1.10.6            
[11] BiocParallel_0.6.1        DESeq2_1.4.5             
[13] RcppArmadillo_0.4.300.8.0 Rcpp_0.11.2              
[15] GenomicRanges_1.16.3      GenomeInfoDb_1.0.2       
[17] IRanges_1.22.8            Biobase_2.24.0           
[19] BiocGenerics_0.10.0      

loaded via a namespace (and not attached):
 [1] annotate_1.42.0    BatchJobs_1.2      BBmisc_1.6         bitops_1.0-6      
 [5] brew_1.0-6         codetools_0.2-8    DBI_0.2-7          digest_0.6.4      
 [9] fail_1.2           foreach_1.4.2      genefilter_1.46.1  geneplotter_1.42.0
[13] grid_3.1.0         hwriter_1.3        iterators_1.0.7    lattice_0.20-29   
[17] locfit_1.5-9.1     plyr_1.8.1         RColorBrewer_1.0-5 RCurl_1.95-4.1    
[21] RSQLite_0.11.4     sendmailR_1.1-2    splines_3.1.0      statmod_1.4.20    
[25] stats4_3.1.0       stringr_0.6.2      survival_2.37-7    tools_3.1.0       
[29] XML_3.98-1.1       xtable_1.7-3       zlibbioc_1.10.0

[Darwin]

I get the exact same error when writing the HTML report. How did you solve this?

> DEXSeqHTML( dxr, FDR=0.1, color=c("#FF000080", "#0000FF80") )
Error in data.frame(..., check.names = FALSE) :
arguments imply differing number of rows: 249794, 0

[Anomilie]

I get the exact same error when writing the HTML report. How did you solve this?

> DEXSeqHTML( dxr, FDR=0.1, color=c("#FF000080", "#0000FF80") )
Error in data.frame(..., check.names = FALSE) :
arguments imply differing number of rows: 249794, 0

I solved this problem by using the transcripts and featureRanges arguments in the DEXSeqDataSet function. For some reason this information doesn't get transferred when using SummarizeOverlaps as the counting method.

The code is as following:

# generate the gtf variable using standard methods
gffFile <- makeTranscriptDbFromGFF("/home/ewilkie/sandrive/bioinformatics/ewilkie/Kenny/Genome_files/Mus_musculus/UCSC/mm10/Annotation/Genes/genes.gtf", format="gtf")
exonicParts <- disjointExons(gffFile , aggregateGenes=FALSE)

## generate transcript annotation
transcripts <- as.list(exonicParts$tx_name)

## generate featureRanges annotation 
t <- as.data.frame(exonicParts)    
test_ranges <- GRanges(seqnames= t$seqnames, IRanges(start=t$start, end=t$end,names=t$gene_id), strand=t$strand)  

## set DEXseqDataSet variable
dxd <- DEXSeqDataSet(raw_data,sampleTable, design, featureID= as.character(exonID), groupID= geneID,  transcripts=transcripts, featureRanges=test_ranges)

## perform the same analysis and call DEXseqHTML as you did before

Hope this helps

[arkanion]

Alternatively, you can specify "flattenedfile" option when reading the dexseq object not to get this error:

DEXSeqDataSetFromHTSeq(countfiles = file.path( inDir, countFiles ), sampleData = sampleData, design = ~ sample + exon + type:exon + condition:exon, [COLOR="Red"]flattenedfile = flattenedfile[/COLOR])

which is the annotation file that was originated with the script "dexseq_prepare_annotation.py".