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 )
DESeq分析流程
by
Tags:
Leave a Reply