SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
[GlimmerHMM] Is my understanding right? syintel87 Bioinformatics 1 02-05-2014 11:10 AM
understanding chromosomal variantions lre1234 General 2 01-08-2014 05:27 AM
Understanding CuffCompare studentofscience RNA Sequencing 0 12-03-2012 10:09 PM
Understanding annotation tahamasoodi Bioinformatics 1 11-19-2012 08:28 AM
help in understanding vinumanikandan General 3 02-10-2011 12:56 AM

Reply
 
Thread Tools
Old 05-06-2014, 09:27 PM   #1
zslee
Member
 
Location: shanghai of china

Join Date: May 2009
Posts: 29
Default dexseq results understanding

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
zslee is offline   Reply With Quote
Old 06-11-2014, 12:25 AM   #2
Anomilie
Member
 
Location: NY

Join Date: Jun 2013
Posts: 13
Default

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?
Anomilie is offline   Reply With Quote
Old 06-17-2014, 05:22 AM   #3
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

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
areyes is offline   Reply With Quote
Old 06-17-2014, 09:35 PM   #4
Anomilie
Member
 
Location: NY

Join Date: Jun 2013
Posts: 13
Default

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.

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")
Anomilie is offline   Reply With Quote
Old 06-20-2014, 06:47 AM   #5
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

you could try dropbox or similar tools!
areyes is offline   Reply With Quote
Old 07-07-2014, 06:06 AM   #6
xuer
Member
 
Location: germany

Join Date: Sep 2008
Posts: 17
Default is it solved?

I have the same problem. I wondering if the problem was solved, if yes, how? could anybody reply?
Thanks!
xuer is offline   Reply With Quote
Old 11-17-2014, 08:27 AM   #7
quinne5
Junior Member
 
Location: dublin

Join Date: Mar 2009
Posts: 5
Default

Just wondering if this has ever been solved -I'm having the same issue..
quinne5 is offline   Reply With Quote
Old 11-18-2014, 06:01 AM   #8
areyes
Senior Member
 
Location: Heidelberg

Join Date: Aug 2010
Posts: 165
Default

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
areyes is offline   Reply With Quote
Old 11-18-2014, 06:52 AM   #9
quinne5
Junior Member
 
Location: dublin

Join Date: Mar 2009
Posts: 5
Default

Thanks for the the update-ill give the development version a go, I have over 20 samples and the analysis is on human genes.
quinne5 is offline   Reply With Quote
Old 09-03-2015, 12:36 AM   #10
wanfahmi
Member
 
Location: North Sea

Join Date: Apr 2008
Posts: 34
Default

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.
wanfahmi is offline   Reply With Quote
Reply

Tags
dexseq

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 11:15 PM.


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