DESeq分析流程

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq")
library("DESeq")
database <- read.table(file = "macrophage_genes.count_table.matrix", sep = "\t", header = T, row.names = 1)
countData <- database[,1:3]
condition <- factor(c("A","B","C"))#无生物学重复
#type <- factor(c(rep("A",3), rep("B",3))) 有生物学重复
database <- round(as.matrix(countData))#取整数型
dds <- newCountDataSet(database,condition)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds, method="blind", sharingMode="fit-only" )#无生物学重复
dds <- estimateDispersions(dds) # 有生物学重复
resAvsB <- nbinomTest(dds,"A","B")
resAvsC <- nbinomTest(dds,"A","C")
table(resAvsB$pval <0.05)
table(resAvsC$pval <0.05)
resAvsB <- resAvsB[order(resAvsB$pval),]#排序
resAvsC <- resAvsC[order(resAvsC$pval),] 
plotMA(resAvsB, ylim = c(-5,5), col = ifelse(resAvsB$pval>=0.05, "gray32", "red3"),linecol = "#ff000080" )#画MA图
norcounts <- counts(dds, normalized=T)#提取标准化后的counts
write.table(as.data.frame(norcounts),
file="normalized_counts_ABC.txt",
sep="\t",
quote = F)
resAvsBup <- subset(resAvsB, pval < 0.05 & log2FoldChange >0 )
resAvsBdown <- subset(resAvsC, pval < 0.05 & log2FoldChange <0 )

发表评论

电子邮件地址不会被公开。 必填项已用*标注