Hi all,
My name is Nandan Deshpande. I am working at the Systems Biology Initiative, UNSW, Sydney.
I have started using DEXSeq recently. I have a few questions about using the tool ...I will appreciate if someone can assist me.
My experiment design (for cancer data-set) includes 9 patients, 2 condition states for each patient ('C1D1' & 'EOT') and 2 batches ('BACTH1' & 'BATCH2') in which the RNASeq data was generated.
Our main question is to check if any of the exons show Differential exon usage across conditions C1D1 and EOT for the nine patients.
I wanted to confirm with the experts in the blog if the steps I have used for my analysis are sufficient to produce a correct result.
My script (which I have tried to derive from the DEXSeq vignette) is as below --->
library("DEXSeq")
library("BiocParallel")
#Using multiple cores
BPPARAM=MulticoreParam(workers=12)flattenedFile <- "XXX_DEXSeq.gff"
sampleTable= data.frame(row.names=c("AW_C1D1","NR_C1D1","DD_C1D1","EL_C1D1","DG_C1D1","EC_C1D1","GR_C1D1","LA_C1D1","RL_C1D1","AW_EOT","NR_EOT","EL_EOT","DD_EOT","EC_EOT","GR_EOT","LA_EOT","RL_EOT","DG_EOT"),condition=c("START","START","START","START","START","START","START","START","START","END","END","END","END","END","END","END","END","END"),libType=c("Batch1","Batch1","Batch2","Batch1","Batch2","Batch2","Batch2","Batch2","Batch2","Batch1","Batch1","Batch1","Batch2","Batch2","Batch2","Batch2","Batch2","Batch2"))
suppressPackageStartupMessages(library("DEXSeq") )
#defining models
formulaFullModel= ~sample+exon+libType:exon+condition:exon
formulaReducedModel= ~sample+exon+libType:exon
dxd=DEXSeqDataSetFromHTSeq(countFile,sampleData=sampleTable,design=~sample+exon+condition:exon,flattenedfile=flattenedFile )
dxd=estimateSizeFactors( dxd )
sizeFactors(dxd)
dxd=estimateDispersions( dxd, formula=formulaFullModel, BPPARAM=BPPARAM )
dxd=testForDEU( dxd,reducedModel= formulaReducedModel,fullModel = formulaFullModel, BPPARAM=BPPARAM)
dxd=estimateExonFoldChanges( dxd,fitExpToVar="condition", BPPARAM=BPPARAM)
dxr2=DEXSeqResults(dxd)
table( dxr2$padj< 0.1)
DEXSeqHTML(dxr2,FDR=0.1,color=c("#FF000080","#0000FF80") )
Other Queries:
- Is there a straight forward way to draw the PCA plot in DEXSeq, using the exon count table? I can then test if I am able to group my patient conditions with and without using the '*libType' in the model (I assume using 'libType' should adjust for variations due to bacth effect!!)
- When using EdgeR/DeSeq for differential expression analysis, we restrict the use of genes for which the cpm values for a certain number of replicates is say >3 cpm; Is there a step similar to this for DEXSeq?
I am getting many exons shown to be having differential exon usage with individual exon counts very small (between 0-10)!! :-(
Regards,
Nandan