Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
library(DESeq2)
#/ from https://raw.githubusercontent.com/shangguandong1996/picture_link/main/WFX_count_Rmatrix.txt
data <- read.delim("WFX_count_Rmatrix.txt", row.names = "Geneid")
get_res <- function(dds){
dds <- DESeq(dds)
res <- results(dds, name=resultsNames(dds)[2])
lfcShrink(dds, coef=resultsNames(dds)[2], res=res, type="ashr")
}
#/ original data:
res1 <-
get_res(
DESeqDataSetFromMatrix(
countData = ceiling(data),
colData = data.frame(condition=factor(rep(c("A", "B"), each=2))),
design = ~condition)
)
#/ prefiltered:
res2 <-
get_res(
DESeqDataSetFromMatrix(
countData = ceiling(data)[rowSums(data > 10 ) >= 2,],
colData = data.frame(condition=factor(rep(c("A", "B"), each=2))),
design = ~condition)
)
#/ the same with simply copy/pasting the samples to make up more replicates,
#/ now with 20 samples per group:
res3 <- get_res(
DESeqDataSetFromMatrix(
countData = ceiling(do.call(cbind, lapply(1:10, function(x) data))),
colData = data.frame(condition=factor(rep(rep(c("A", "B"), each=2), 10))),
design = ~condition)
)
#/ MA-plots are essentially the same towards odd shape,
#/ color is turned off here to not look messy
pdf("maplots.pdf")
par(mfrow=c(2,2))
plotMA(res1, alpha=0, ylim=c(-10,10), main="original")
plotMA(res2, alpha=0, ylim=c(-10,10), main="prefiltered")
plotMA(res3, alpha=0, ylim=c(-10,10), main="20 samples per group")
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment