Entering edit mode
Hello Everyone,
I am in a big problem, when I am running the DESeq2, the genes which are up and down-regulated or at all differentially expressed, some of them are not appearing in my result when I am adjusting the p-value <0.05, and when I am checking before adjusting they appear because their pvalue are more than that like around 0.3 or 0.1, I need your advice and suggestion to find the problem? also, I am sharing the commands which I use for running the DESeq2. please let me know what is the problems? I would really appreciate your help in my this tough situation.
require(DESeq2)
Counts=read.table("~/Desktop/count/SRR00000.count")[,2]
for (i in dir("~/Desktop/count/", full.names = T)[-1]){
Counts=cbind(Counts,read.table(i)[,2])
}
rownames(Counts)=read.table("~/Desktop/count/SRR00000.count")[,1]
colnames(Counts)=gsub(".count","",dir("~/Desktop/count/"))
Counts=Counts[-grep("__",rownames(Counts)),]
countsMatrix=as.matrix(Counts)
class=read.csv("~/Desktop/Meta-analysis/table.csv")[,2]
class=data.frame(condition=as.factor(class))
for (i in 1:ncol(Counts)){
Counts[,i]=Counts[,i]/CountsNorm$sizeFactor[i]
}
carrier=rowMeans(Counts[,which(CountsNorm$condition=="carrier")])
noncarrier=rowMeans(Counts[,which(CountsNorm$condition=="noncarrier")])
FoldChange=carrier/noncarrier
log2FoldChange=log2(FoldChange)
Pvalue=c()
carrierindex=which(CountsNorm$condition=="carrier")
noncarrierindex=which(CountsNorm$condition=="noncarrier")
for (i in 1:nrow(Counts)){
Pvalue=c(Pvalue,t.test(Counts[i,carrierindex],Counts[i,noncarrierindex])$p.value)
}
UpRegulateGenes=which(log2FoldChange>=1 & Pvalue<0.05)
DownRegulatedGenes=which(log2FoldChange<=-1 & Pvalue<0.05)
results=data.frame(log2FoldChange, Pvalue)
DEG=which(log2FoldChange>=1 & Pvalue<0.05 | log2FoldChange<=-1 & Pvalue<0.05)
results=data.frame(log2FoldChange,Pvalue)
results=results[order(results$Pvalue),]
results=results[which(results$Pvalue<0.05),]
results=results[order(results$log2FoldChange),]
write.table(results, file = "result.txt", row.names = T, col.names = T, sep = "\t")
Sorry, not got the question properly?
You did not run DESeq2 for DE.
you mean this?
results <- DESeq(dds)
Connecting this thread to the new thread:
t.test to determine the P-value in DESeq2
Many thanks for your advise, I tried these commands for running the DESeq2: ,,,
library(DESeq2)
count_matrix <- as.matrix(read.csv("~/Desktop/sample.csv", row.names = 1))
head(count_matrix, 2)
View(count_matrix)
condition <- factor(c(rep("carrier",3), rep("noncarrier", 3)))
sampleTable <- data.frame(condition = as.factor(condition))
rownames(sampleTable) <- colnames(count_matrix)
sampleTable
deseq <- DESeqDataSetFromMatrix(countData = count_matrix,
deseq$condition <- relevel(dds$condition, ref = "noncarrier")
d.deseq <- DESeq(deseq)
res<- results(d.deseq)
res<- results(d.deseq, alpha = 0.05 )
summary(res)
- but the number of DEGs are very less as below:
out of 32159 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up) : 23, 0.072%
LFC < 0 (down) : 10, 0.031%
outliers [1] : 168, 0.52%
low counts [2] : 2542, 7.9%
(mean count < 0)
it gives me these below command too,
[1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results
what do you think where the problem might be? I would really appreciate your help to run DESeq2 successfully, I tried many commands but still the same problem.
The first thing to be 100% sure of is that you've got your conditions applied to the correct samples, but if that's right, then I don't see why you would not believe these results.
Many thanks for your answer. Because I am recreating the data from the original paper. In that paper, the number of genes that are up and down-regulated is higher than my results. That's why I do not believe this results. what do you recommend for me in this situation? as I read some manuals it was written that for less than 12 replicates it is better to use edge R, and for more than 12 using DESeq2 is better. Do you think it is because of that I can not get the DEGs as of the original papers?