DESeq2和TopGO的R代码

library(DESeq2)

library(ggplot2)

#library(ggrepel)

raw.data1<-read.table("feature_count_combo4DESeq2.tab", check.names = FALSE, header=T, row.names = 1)

#raw.data1 <- raw.data1[ rowSums(raw.data1) > 1, ]

#dim(raw.data1)

#raw.data1 <- replace(raw.data1, raw.data1 == 0, 1)

idx_col = grep("FN_45|LN_45", colnames(raw.data1), perl=TRUE)

raw.data_sub <- raw.data1[,idx_col]

sam_info <- gsub("_[12]$", "", colnames(raw.data_sub), perl = T)

colData  <- data.frame(LN = sam_info)

rownames(colData) = colnames(raw.data_sub)

dds <- DESeqDataSetFromMatrix(countData = raw.data_sub,

                              colData = colData,

                              design =~ LN)

colData(dds)$LN = factor(colData(dds)$LN,levels=c("FN_45", "LN_45"))

dds <- DESeq(dds)

res <- results(dds)

## Generate the table

resOrd<-res[order(res$padj),]

res1_table<-as.data.frame(resOrd)

res1_table$fc<-2^res1_table$log2FoldChange ##generate the FC column

res1_table$gene_id = rownames(res1_table)

res1_table<-res1_table[,c(8,1,2,7,3,4,5,6)]

res1_table_noNA <- res1_table[complete.cases(res1_table), ]

res1_table_sig <- res1_table_noNA[res1_table_noNA$padj < 0.05, ]

geneid_sig <- gsub("-RA$", "", rownames(res1_table_sig), perl=T)

library(topGO)

geneID2GO <- readMappings("/media/xie186/easystore/WIC/project/danshen_RNA_seq_metabolomcis/data/GO_KEGG_info/salvia_cds_GO4topGO.tab")

GO2geneID <- inverseList(geneID2GO)

gene_number_GO = do.call(rbind, lapply(GO2geneID, function(x) length(x)))

geneNames <- names(geneID2GO)

geneList <- factor(as.integer(geneNames %in% geneid_sig))

names(geneList) <- geneNames

GOdata <- new("topGOdata", ontology = "MF", allGenes = geneList,

              nodeSize = 5,  #It is often the case that many GO terms which have few annotated genes are detected to be significantly

              #enriched due to artifacts in the statistical test. These small sized GO terms are of less importance for the

              #analysis and in many cases they can be omitted. By using the nodeSize argument the user can control the

              #size of the GO terms used in the analysis. Once the genes are annotated to the each GO term and the true

              #path rule is applied the nodes with less than nodeSize annotated genes are removed from the GO hierarchy.

              #We found that values between 5 and 10 for the nodeSize parameter yield more stable results. The default

              #value for the nodeSize parameter is 1, meaning that no pruning is performed.

              annot = annFUN.gene2GO, gene2GO = geneID2GO)

### The list of genes of interest can be accessed using the method sigGenes():

#sg <- sigGenes(GOdata)

#str(sg)

#numSigGenes(GOdata)

resultFisher <- runTest(GOdata, algorithm = "classic", statistic = "fisher")

resultKS <- runTest(GOdata, algorithm = "classic", statistic = "ks")

resultKS.elim <- runTest(GOdata, algorithm = "elim", statistic = "ks")

#allRes <- GenTable(GOdata, classicFisher = resultFisher,

#                    classicKS = resultKS, elimKS = resultKS.elim,

#                    orderBy = "elimKS", ranksOf = "classicFisher", topNodes = 10)

allRes <- GenTable(GOdata,

                  classicFisher = resultFisher,

                  #classicKS = resultKS,

                  #elimKS = resultKS.elim,

                  orderBy = "classicFisher", ranksOf = "elimKS", topNodes = length(attributes(resultKS.elim)$score))

#We can visualise the position of the statistically significant GO terms

#in the GO hierarchy by using the following functions:

showSigOfNodes(GOdata, score(resultKS.elim), firstSigNodes = 5, useInfo = 'all')

#The second command makes a pdf file ("tGO_classic_5_all.pdf") with the picture.

#The significant GO terms are shown as rectangles in the picture. The most significant

#terms are coloured red and least significant in yellow:

printGraph(myGOdata, resultFisher, firstSigNodes = 5, fn.prefix = "tGO", useInfo = "all", pdfSW = TRUE)

myterms <- allRes$GO.ID

mygenes <- genesInTerm(GOdata, myterms)

allRes$gene.list = rep(NA, length(allRes$GO.ID))

for (i in 1:length(myterms))

{

  myterm <- myterms[i]                      ## get the term

  mygenesforterm <- mygenes[myterm][[1]]    # get the genes list

  myfactor <- mygenesforterm %in% geneid_sig # find the genes that are in the list of genes of interest

  mygenesforterm2 <- mygenesforterm[myfactor == TRUE]    ## get the genes

  mygenesforterm2 <- paste(mygenesforterm2, collapse=',')

  allRes$gene.list[i] = mygenesforterm2

  #print(paste("Term",myterm,"genes:",mygenesforterm2))

}