I'm clustering genes in terms of expression after clustering them. I'm taking out clusters and trying to find out what kind of GO terms are coming up. I came across this figure from this paper, I want to do as the figure B:
Here you can find the file I intend to use to make something like the above plot.
I am using this the code from a question in Biostar:
DAVIDfile <- "DAVID.txt"
DAVID <- read.table(DAVIDfile, sep="\t", header=TRUE)
colnames(DAVID)
class(DAVID)
enrichBcutoff <- 0.05
DAVID <- subset(DAVID, Benjamini< enrichBcutoff)
DAVID <- DAVID[,c(1,2,3,5,6,10)]
head(DAVID)
A <- data.frame(sigGeneList)
annGSEA <- data.frame(row.names=A)
for (j in 1:length(sigGeneList))
{
gene <- sigGeneList[j]
pattern <- paste("^", gene, ", |, ", gene, "$| ", gene, ",", sep="")
for (k in 1:nrow(DAVID))
{
if (any(grepl(pattern, DAVID$Genes[k])))
{
annGSEA[j,k] <- 1
}
else
{
annGSEA[j,k] <- 0
}
}
}
colnames(annGSEA) <- DAVID[,2]
annGSEA <- annGSEA[,apply(annGSEA, 2, mean)!=0]
annGSEA <- annGSEA[apply(annGSEA, 1, mean)!=0,]
require(ComplexHeatmap)
require(circlize)
topTable <- topTable[which(topTable$Gene %in% rownames(annGSEA)),]
topTable <- topTable[match(rownames(annGSEA), topTable$Gene),]
geneLab=10
termLab=8
#Create heatmap annotations
dfMinusLog10FDRGenes <- data.frame(-log10(topTable[which(topTable[,1] %in% rownames(annGSEA)),"padj"]))
dfMinusLog10FDRGenes[dfMinusLog10FDRGenes=="Inf"] <- 0
dfFoldChangeGenes <- data.frame(topTable[which(topTable[,1] %in% rownames(annGSEA)),"log2FoldChange"])
dfGeneAnno <- data.frame(dfMinusLog10FDRGenes, dfFoldChangeGenes)
colnames(dfGeneAnno) <- c("DEG\nsignificance\nscore", "Regulation")
dfGeneAnno[,2] <- ifelse(dfGeneAnno[,2]>0, "Up-regulated", "Down-regulated")
colours <- list("Regulation"=c("Up-regulated"="royalblue", "Down-regulated"="yellow"))
haGenes <- rowAnnotation(df=dfGeneAnno, col=colours, width=unit(1,"cm"))
dfMinusLog10BenjaminiTerms <- data.frame(-log10(read.table(DAVIDfile, sep="\t", header=TRUE)[which(read.table(DAVIDfile, sep="\t", header=TRUE)$Term %in% colnames(annGSEA)),"Benjamini"]))
colnames(dfMinusLog10BenjaminiTerms) <- "GO Term\nsignificance\nscore"
haTerms <- HeatmapAnnotation(df=dfMinusLog10BenjaminiTerms,
colname=anno_text(colnames(annGSEA), rot=40, just="right", offset=unit(1,"npc")-unit(2,"mm"), gp=gpar(fontsize=termLab)),
annotation_height=unit.c(unit(1, "cm"), unit(8, "cm")))
pdf("GO.pdf", width=7, height=12)
hmapGSEA <- Heatmap(annGSEA,
name="My enrichment",
split=dfGeneAnno[,2],
col=c("0"="white", "1"="forestgreen"),
rect_gp=gpar(col="grey85"),
cluster_rows=T,
show_row_dend=T,
row_title="Statistically-significant genes",
row_title_side="left",
row_title_gp=gpar(fontsize=12, fontface="bold"),
row_title_rot=0,
show_row_names=TRUE,
row_names_gp=gpar(fontsize=geneLab, fontface="bold"),
row_names_side="left",
row_names_max_width=unit(15, "cm"),
row_dend_width=unit(10,"mm"),
cluster_columns=T,
show_column_dend=T,
column_title="Enriched terms",
column_title_side="top",
column_title_gp=gpar(fontsize=12, fontface="bold"),
column_title_rot=0,
show_column_names=FALSE,
#column_names_gp=gpar(fontsize=termLab, fontface="bold"),
#column_names_max_height=unit(15, "cm"),
show_heatmap_legend=FALSE,
#width=unit(12.5, "cm"),
clustering_distance_columns="euclidean",
clustering_method_columns="ward.D2",
clustering_distance_rows="euclidean",
clustering_method_rows="ward.D2",
bottom_annotation=haTerms)
draw(hmapGSEA + haGenes, heatmap_legend_side="right", annotation_legend_side="right")
dev.off()
However there is an error message after annGSEA:
dput(annGSEA)
structure(list(), .Names = character(0), row.names = "c(\"LRP1\", \"ELAVL4\", \"ENOX1\", \"TRIP6\", \"KHDRBS3\", \"TDRD9\", \"THOC6\", \"CPSF1\", \"HSPB1\", \"SF3A2\", \"STXBP1\", \"DUS3L\", \"INTS1\", \"DZIP1L\", \"HIST1H1B\", \"TRMT2A\", \"RPUSD1\", \"MRPS2\", \"EXOSC5\", \"TDRD1\", \"PURG\", \"NOL12\", \"TRMT61A\", \"XAB2\", \"TRMT1\", \"PTRF\", \"POLR2I\", \"MAZ\", \"TBL3\", \"BOP1\", \"RBM10\", \"TSR3\", \"SART1\", \"CCDC86\", \"SAMD4A\", \"RBM11\", \"RPS5\", \"HNRNPLL\", \"POLRMT\", \"NSUN7\", \"PABPC5\", \"DDX54\", \"PRRC2B\", \"RPS15\", \"MYBBP1A\", \"DUS1L\", \"ARHGEF28\", \"RPL28\", \"ZNF598\", \"ZMAT5\", \"RPL36\", \"RPL8\", \"DHX37\", \"ZC3H7B\", \n\"RRP9\", \"PQBP1\", \"ZNF385A\", \"PUS1\", \"DDX49\", \"TCOF1\", \"BICC1\", \"GLTSCR2\", \"SCAF1\", \"MRPL4\", \"FSCN1\", \"EDF1\", \"NOVA1\", \"PCBP3\", \"FBL\", \"MRPL55\", \"PUF60\", \"NYNRIN\", \"ZFP36L1\", \"RPS9\", \"MKRN3\", \"GTF2F1\", \"EEFSEC\", \"FARSA\", \"TARBP2\", \"RPL35\", \"RPS19\", \"ZC3H3\", \"H1F0\", \"BST2\", \"FASN\", \"ZYX\", \"NOC4L\", \"LONP1\", \"RUVBL2\", \"MRPL38\", \"TRAP1\", \"TUFM\", \"RBM42\", \"TNS1\", \"EDC4\", \"WDR6\", \"PCBP4\", \"PUSL1\", \"EEF2\", \"RPS28\", \"RAVER1\", \"APEH\", \"TSEN54\", \"VARS\", \"C7orf50\", \"GRN\", \"ARHGEF1\", \"SF3B5\", \"RPL18\", \"EIF3G\", \n\"ANG\", \"EIF5AL1\", \"HIST1H1C\", \"RPL13\", \"U2AF1L4\", \"C9orf114\", \"FASTK\", \"TSEN34\", \"SRRT\", \"RPS2\", \"MRPL54\", \"CSRP1\", \"HELZ2\", \"RPLP2\", \"RPL18A\", \"PRRC2A\", \"GTPBP1\", \"NOC2L\", \"RBM19\", \"RNASEH2A\", \"ADARB2\", \"SERPINH1\", \"CCDC9\", \"DXO\", \"GNB2L1\", \"FLNA\", \"RNASEH2C\", \"IPO13\", \"HABP4\", \"C11orf68\", \"QTRT1\", \"DGCR14\", \"RPL3\", \"SETD1A\", \"GEMIN4\", \"NPM2\", \"MRPS26\", \"RRP7A\", \"CCDC124\", \"CACTIN\", \"ISG20\", \"INTS5\", \"SNTB2\", \"RPUSD3\", \"RPL29\", \"PPRC1\", \"PPAN\", \"PRPF40B\", \"RBFOX2\", \"C15orf52\", \"ALYREF\", \"EXOSC4\", \n\"GNL1\", \"SIDT2\", \"REPIN1\", \"CHERP\", \"NOB1\", \"RPLP1\", \"EIF5A2\", \"RRP1\", \"MRPL41\", \"RPL19\", \"SUGP1\", \"NAA38\", \"CDC42EP4\", \"EEF1D\", \"DDX41\", \"PRPF19\", \"CLASRP\", \"PRPF6\", \"ZC3H12B\", \"SNRPB\", \"MEPCE\", \"COL14A1\", \"WBSCR16\", \"LSM7\", \"PELP1\", \"VARS2\", \"TRPT1\", \"RPS10\", \"PTBP1\", \"RPS8\", \"CPSF3L\", \"NOSIP\", \"TCF20\", \"PTCD1\", \"APEX1\", \"C1orf35\", \"MRPL23\", \"TWF2\", \"CTIF\", \"DHX38\", \"MRPS34\", \"SYMPK\", \"U2AF2\", \"QARS\", \"CCAR2\", \"RPS3\", \"RPS16\", \"P4HB\", \"MVP\", \"POLR2J\", \"SKIV2L\", \"GTPBP3\", \"ZNF768\", \"UBA52\", \"SUPT5H\", \n\"POLR2L\", \"POLR2A\", \"EIF3K\", \"R3HCC1\", \"PIH1D1\", \"ATXN2L\", \"OAS3\", \"IPO4\", \"CCDC137\", \"SETD1B\", \"ILF3\", \"TBRG4\", \"PKM\", \"AARS2\", \"FAU\", \"ACTN4\", \"RPS19BP1\", \"ANGEL1\", \"PLD6\", \"MOV10\", \"EIF6\", \"DRG2\", \"AARS\", \"PRR3\", \"AC004381.6\", \"RBPMS\", \"CPSF4\", \"NOL6\", \"RPL13A\", \"NPM3\", \"UBA1\", \"REXO4\", \"RRP12\", \"UPF1\", \"RPL11\", \"DDX43\", \"SNRPA\", \"ISY1\", \"MRPS18A\", \"RBMS1\", \"CSTB\", \"FKBP3\", \"MTIF3\", \"MRPL39\", \"CLK4\", \"SRP54\", \"ERI1\", \"CLK1\", \"IMP3\", \"GRSF1\", \"GTF3A\", \"FASTKD1\", \"DYNC1LI1\", \"PTGES3\", \"TRIT1\", \"EXOSC3\", \n\"IGF2BP3\", \"TRIM71\", \"MCTS1\", \"GEMIN2\", \"SUMO1\", \"NQO1\", \"HSPA8\", \"POLR2G\", \"TAF9\", \"MRPL33\", \"FAM46A\", \"MRPS36\", \"FNDC3B\", \"DDX52\", \"DYNLL1\", \"RBM43\", \"MANF\", \"TDRD7\", \"AHNAK\", \"IFIT2\", \"ZMAT3\", \"APOBEC3H\", \"RTN4\", \"S100A4\", \"PIH1D2\", \"ANXA2\", \"PIWIL3\", \"TMSB4X\", \"PIWIL4\", \"RNASE6\", \"LGALS1\", \"RNASE3\", \"TLR7\", \"MAEL\", \"SAMHD1\", \"RNASE2\", \"C2orf15\", \"SIDT1\", \"HERC5\", \"PSIP1\", \"RPL17\", \"ZC3H12C\", \"SLC16A3\", \"TEP1\", \"HEXIM1\", \"CPEB4\", \"ZC3H12D\", \"CORO1A\", \"ZFP36\", \"KIF1C\", \"WARS\", \"CCRN4L\", \"AFF2\", \n\"TST\", \"RBM47\", \"LGALS3\", \"NXF3\", \"PPARGC1A\", \"TLR8\")", class = "data.frame")
>
pattern
, the error points to there being something wrong with that. Also, it looks like you didn't copy paste the error message correctly. Give us the full message and it'll be easier the pin the issue down. $\endgroup$