Hi. I am a 'R' starter. And recently, I analysed GEO database about my research fields. I have a question about results of my analysis. For example, I analysed GEO data which numbered GSE92506. This data shows about the reaction of HUVEC cell to disturbed shear stress model. They opened some analytic data about their experiment and showed me expression profiling by high throughput sequencing. Their data analysis were conducted with HTseq and DESeq2. I downloaded their RNAseq data and analysed with 'R' using DESeq2 protocol and check the DEG level by variable(laminar shear stress(abbreviated "PS"), disterbed shear stress(abbreviated "OS"). Most of my analytic data were similar with their data but fold change level and p value differed slightly. And in case of some gene which are not differed (low level of fold change), my and their data were differed much.
Please tell me about my mistake on my 'R' code. Please~~ㅜㅜ
#DEseq analysis
library(DESeq2)
ddsFullCountTable <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ flow.type)
dds <- ddsFullCountTable[ , ddsFullCountTable$flow.duration == "24h" ]
dds$flow.duration <- droplevels( dds$flow.duration )
dds$flow.type <- relevel( dds$flow.type, "PS" )
dds <- DESeq(dds)
res <- results( dds )
rld <- rlog( ddsFullCountTable )
degcountdata = as.data.frame(assay(rld))
#DEG
eset=ExpressionSet(assayData = as.matrix(degcountdata), featureData = AnnotatedDataFrame(featureData), phenoData = AnnotatedDataFrame(coldata))
fData(eset)=data.frame(fData(eset), res$log2FoldChange, res$pvalue)
deg.idx=which(abs((fData(eset)$res.log2FoldChange)) > 2.0 & (fData(eset)$res.pvalue) < 0.01)
deg.eset=eset[deg.idx,]
as.data.frame(deg.eset)
deg.esetdata = as.data.frame(deg.eset)
fData(deg.eset)
fdat=fData(deg.eset)
row.od=order((fData(deg.eset)$res.log2FoldChange))
fdat=(fData(deg.eset))[row.od,]
write.csv( as.data.frame(fdat), file="f24dat.csv" )
I look forward to your answers. Thank you for reading long question.
