Question: edgeR steps for gene expressions
1
gravatar for Tania
2.4 years ago by
Tania120
Tania120 wrote:

Hello Everyone

I am new to edgeR. I am expecting weird gene expressions. I want to double check my code steps before discovering where else the problem could be. Do these steps make sense to you?

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,ignoreTxVersion = TRUE)
write.csv(txi$counts, "counts.csv")
cts <- txi$counts

First 10 samples are control second 10 samples are tumor:

group = factor(c(rep("Control", 10), rep("Tumor",10)) ) 
dge = DGEList(counts=cts, genes= rownames(data), group=group)
countsPerMillion <- cpm(dge)
summary(countsPerMillion)
countCheck <- countsPerMillion > 1
summary(countCheck)
keep <- which(rowSums(countCheck) >= 2)
dge <- dge[keep,]
summary(cpm(dge))
dge <- calcNormFactors(dge, method="TMM")
dge <- estimateCommonDisp(dge)
dge <- estimateTagwiseDisp(dge)
dge <- estimateTrendedDisp(dge)
et <- exactTest(dge, pair=c("Control", "Tumor"))
etp <- topTags(et, n=100000)
write.csv(etp$table, "dge.csv")

Thanks

rna-seq • 957 views
ADD COMMENTlink modified 2.4 years ago by theobroma221.1k • written 2.4 years ago by Tania120

I am expecting weird gene expressions.

So are you getting normal results instead :)

On a serious note, what is unexpected about the results you are getting?

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by genomax80k

=D Sorry I mean found weird results :) what I expected to be enriched in tumor are enriched in control. I want to make sure code is right before investigating more.

ADD REPLYlink written 2.4 years ago by Tania120

Why didn't you just use estimateDisp? It estimates Common, tagwise and trended dispersions in one run :D, also I would change the value of topTags, since the default value seems to be 1, you would really be getting a list of the first 100,000 genes sorted by pvalue.

ADD REPLYlink written 2.4 years ago by biofalconch440

Ok, Thanks biofalconch :) Is this line correct?

"group = factor(c(rep("Control", 10), rep("Tumor",10)) ) "
ADD REPLYlink written 2.4 years ago by Tania120
1

It seems alright to me, it really depends on the counts matrix though. A PCA might be useful in this case

ADD REPLYlink written 2.4 years ago by biofalconch440
0
gravatar for theobroma22
2.4 years ago by
theobroma221.1k
theobroma221.1k wrote:

Something doesn’t seem right to me up top...did you paste all of your code?

I could be wrong but where is your “data” object for the “dge” object in:
dge = DGEList(counts=cts, genes= rownames(data), group=group)...?

ADD COMMENTlink written 2.4 years ago by theobroma221.1k

its in cts in the above part from tximport? Wrong?

ADD REPLYlink written 2.4 years ago by Tania120
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1052 users visited in the last hour