dexseq results understanding

[zslee]

Hey, Now I am trying to find splicing differences between two tumor and two normal samples(so only one condition: N and T factors), and prepare data following your vignette, everything goes well, but after the last step(bellow), I got the dxr dataframe I find amost half lines with exon usage coefficient as NA, when I extract the corresponding gene I find the exons well covered (one example bellow), so, my question is why so many NAs? could anyone please give some possible explanations which I would check.
Many Many Thanks.

######
my processes:
dxd = DEXSeqDataSetFromHTSeq(countFiles, sampleData=sampleTable,design=~sample + exon + condition:exon, flattenedfile=flattenedFile)
dxd=estimateSizeFactors(dxd)
dxd=estimateDispersions(dxd)
dxd=testForDEU(dxd)
dxd=estimateExonFoldChanges(dxd, fitExpToVar="condition")
dxr=DEXSeqResults(dxd)

example:

head(dxr)

LRT p-value: full vs reduced

DataFrame with 6 rows and 13 columns
groupID featureID exonBaseMean dispersion stat pvalue padj R S log2fold_R_S genomicData countData transcripts
<character> <character> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <GRanges> <matrix> <list>
ENSG00000000003:E001 ENSG00000000003 E001 259.50 0.005893988 0.7015350362 0.4022684 1 NA NA NA chrX:-:[99883667, 99884983] 199 201 471 ... ########
ENSG00000000003:E002 ENSG00000000003 E002 112.00 0.001052854 1.7011996938 0.1921312 1 NA NA NA chrX:-:[99885756, 99885863] 98 96 183 ... ########
ENSG00000000003:E003 ENSG00000000003 E003 92.75 0.005572201 0.0007320773 0.9784143 1 NA NA NA chrX:-:[99887482, 99887537] 97 76 145 ... ########
ENSG00000000003:E004 ENSG00000000003 E004 71.50 0.010739013 0.0385784467 0.8442862 1 NA NA NA chrX:-:[99887538, 99887565] 76 64 112 ... ########
ENSG00000000003:E005 ENSG00000000003 E005 81.00 0.007045578 0.7014820356 0.4022862 1 NA NA NA chrX:-:[99888402, 99888438] 90 61 131 ... ########
ENSG00000000003:E006 ENSG00000000003 E006 111.75 0.002231303 1.9566454114 0.1618725 1 NA NA NA chrX:-:[99888439, 99888536] 117 80 192 ... ########

exon counts:
head(counts(dxd))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
ENSG00000000003:E001 199 201 471 167 855 737 1587 511
ENSG00000000003:E002 98 96 183 71 956 842 1875 607
ENSG00000000003:E003 97 76 145 53 957 862 1913 625
ENSG00000000003:E004 76 64 112 34 978 874 1946 644
ENSG00000000003:E005 90 61 131 42 964 877 1927 636
ENSG00000000003:E006 117 80 192 58 937 858 1866 620

[Anomilie]

I'm having a similar issue where the log2fold value for many exons is also NA, while they show sufficient counts in the data. Did you find a solution/explanation for this behavior?

[areyes]

Thanks for reporting this! The current implementation of DEXSeq estimates fold changes only for those exons for which a p-adjustes value was calculated (e.g. not being an outlier).

How many of such cases of having a p-value but not a fold change do you see in your data? If they are many, could you maybe send me your object so I could have a closer look to what is happening?

Alejandro

[Anomilie]

Thanks for your response.

Out of the total 716 genes that are significant (FDR < 0.05), 101 have NA for the log2fold column, but they all have an adjusted p-value.

I am not able to attach the .Rdata object of the dxd variable created using the documentation as the file is too large. Can I send it to you in any other way or is there a more specific object that would be useful for you to debug what is going on? I generated the dxd object using the following code.

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

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


align <- "/Samples"

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

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

genehits <- summarizeOverlaps(exonicParts, bf1, mode="IntersectionStrict", ignore.strand=FALSE, singleEnd=FALSE, inter.feature=TRUE, 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


dxd <- DEXSeqDataSet(raw_data,sampleTable, design, featureID= as.character(exonID), groupID= g)

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

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

save(dxd, file="DEXseq_v1.Rdata")

[areyes]

you could try dropbox or similar tools!

[xuer]

is it solved?

I have the same problem. I wondering if the problem was solved, if yes, how? could anybody reply?
Thanks!

[quinne5]

Just wondering if this has ever been solved -I'm having the same issue..

[areyes]

Sorry for the very late reply. Could you confirm if this happens when you have a moderate number of samples (around 10-12) and for genes with lots of exonic regions?

I think this has to do with the GLMs that are calculated for each gene to estimate the exon fold changes. The way this works is that a model frame is created where the rows is number of exonic regions times number of samples. I set up a threshold in 3000 to the number of rows of the model frame such the fit would not be done if the model frame for a gene passes this threshold. I added an option in the latest development version such that users can increase this number, but consider that the larger this value, the larger it will take to compute. The parameter is the maxRowsMF of the function estimateExonFoldChanges.

Consider that this is a temporary solution, we need to think in a smarter way to deal with this!

Alejandro

[quinne5]

Thanks for the the update-ill give the development version a go, I have over 20 samples and the analysis is on human genes.

[wanfahmi]

Can DEXSeq handle more than 70 samples to analyse? I have different human tissue to analyse which not same as normal vs disease as shown in vignette.