转录本差异DRIMSeq使用

官方文件
参考1
参考2

rm(list = ls())
suppressMessages(library(DRIMSeq))
counts<-read.csv("/your_dir/tables/count.cluster.csv",header = T,row.names = 1)
new_counts<- data.frame(gene_id=unlist(lapply(row.names(counts), FUN = function(x) {return(strsplit(x, split = "_",fixed = T)[[1]][1])})),
                                 feature_id=row.names(counts),
                                 counts)

samples<-data.frame(sample_id=colnames(counts),
                    condition=c(rep("fed",3),rep("fasting",3)))
d <- dmDSdata(counts=new_counts, samples=samples)
#methods(class=class(d))
# filter the object, adjust the parameters
n <- 6 # total samples
n.small <- 3 # the samples numbers of the smallest group
d_pass <- dmFilter(d,
              min_samps_feature_expr=n.small, min_feature_expr=10,
              min_samps_feature_prop=n.small, min_feature_prop=0.1,
              min_samps_gene_expr=n, min_gene_expr=10)


design_full <- model.matrix(~condition, data=DRIMSeq::samples(d_pass))
colnames(design_full)
d <- dmPrecision(d_pass, design=design_full)
d <- dmFit(d, design=design_full)
d <- dmTest(d, coef="conditionfed")

res <- DRIMSeq::results(d)
head(res)

res.txp <- DRIMSeq::results(d, level="feature")
head(res.txp)


no.na <- function(x) ifelse(is.na(x), 1, x)
res$pvalue <- no.na(res$pvalue)
res.txp$pvalue <- no.na(res.txp$pvalue)

res<-res[order(res$adj_pvalue),]
res.txp<-res[order(res.txp$adj_pvalue),]
idx <- which(res$adj_pvalue < 0.05)[1]
res[idx,]

plotProportions(d, res$gene_id[idx], "condition")
save.image(file = "/your_dir/DRIM_results/DRIM_seq.RData")
write.csv(res, "/your_dir/DRIM_results/gene_DRIM_Seq.csv")
write.csv(res.txp, "/your_dir/DRIM_results/feature_DRIM_Seq.csv")