Was this page helpful?

DESeq2

    Table of contents
    No headers
    ##load the library
    library(DESeq2)
     
    ##Set the path where the data file is
    setwd("C:/Users/Test/")
     
    ##read in the data file
    raw.data <- read.delim("Final_MergedCounts.txt", header=TRUE,stringsAsFactors=FALSE)
    head(raw.data)
    dim(raw.data)
    raw.data<-na.omit(raw.data)
     
    ##replace all the "0" to "1"
    raw.data <- replace(raw.data, raw.data == 0, 1)
    head(raw.data)
     
    ##input all the count data
    countData <- as.matrix(raw.data[,2:13],genotype = factor(c(
                        "G1","G1","G1",    "G1","G1","G1","G2","G2","G2","G2","G2","G2")),
     
                          time=factor(c("T0","T0","T0", "T1","T1","T1",T0","T0","T0", "T1","T1","T1")))
    head(countData)
     
    ##input the factor levels
    ExpDesign <- data.frame(row.names=colnames(raw.data)[2:61],
      genotype=factor(c("G1","G1","G1",    "G1","G1","G1","G2","G2","G2","G2","G2","G2")),
     
      time=factor(c("T0","T0","T0", "T1","T1","T1", 
                    "T0","T0","T0", "T1","T1","T1",)))
     
    head(ExpDesign)
     
     
     
    ##############################################################################################
     
    #######  Use design matrix with two main factors and the interaction######################
     
    ##############################################################################################
     
    CDS <- DESeqDataSetFromMatrix(countData,colData=ExpDesign,design=~genotype+time+genotype:time)
    rownames(CDS)<- raw.data$ID
    CDS
     
    ##Test results by contrasting
    CDS <- DESeq(CDS)
     
    resultsNames(CDS)
     
     
     
    ##get the results for the contrasts of interest and put into csv file
    ##Compare within each genotype
    ##note:the first element is -1 and the second element is +1
    res <- results(CDS,contrast=list("genotypeG1.timeT1","genotypeG1.timeT0"))
    head(res,5)
    write.csv(as.data.frame(res), file="C:/Users/result.csv")
    setwd("C:/Users/")
    ##MA-plot
    ##Need to order the results by padj for the significant genes (padj<0.1) to be shown as red
    resOrdered<-res[order(res$padj),]
    head(resOrdered)
    plotMA(res,main="DESeq2", ylim=c(-2,2))
     
     
    ##Compare within each time point, other genotype vs B73
    ##get the results for the contrasts of interest and put into csv file
    ##Compare within each genotype
    ##note:the first element is -1 and the second element is +1
    res <- results(CDS,contrast=list("genotypeG1.timeT0","genotypeG2.timeT0"))
    head(res,5)
    write.csv(as.data.frame(res), file="C:/Users/result2.csv")
    ##MA-plot
    ##Need to order the results by padj for the significant genes (padj<0.1) to be shown as red
    resOrdered<-res[order(res$padj),]
    head(resOrdered)
    setwd("C:/Users/")
    plotMA(res,main="DESeq2", ylim=c(-2,2))
     
     
     
    Was this page helpful?
    Tag page (Edit tags)
    • No tags
    You must login to post a comment.