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))
}
