2
$\begingroup$

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:the second figure

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")
>
Improve this question
$\endgroup$
6
  • $\begingroup$ I think you forgot the link to the paper you are referencing. You tagged the question with r, which problem have with which code? Perhaps you want to make the figure B and you don't know how? $\endgroup$
    – llrs
    Oct 31, 2018 at 9:53
  • $\begingroup$ yes i m trying to make figure and truly i don;t know how to begin with $\endgroup$
    – kcm
    Oct 31, 2018 at 10:20
  • $\begingroup$ You added a lot of code, you already seem how to start! What is the output of this code? (could you upload the hmapGSEA plot?) What didn't you like or is not working as expected? $\endgroup$
    – llrs
    Oct 31, 2018 at 10:31
  • $\begingroup$ I took this code from biostar, im stuck at the annGSEA after that when i loop through it doesn't match the pattern of genes would give the error message in the question have a look $\endgroup$
    – kcm
    Oct 31, 2018 at 13:13
  • $\begingroup$ Please provide the output for the variable 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$
    – Kohl Kinning
    Oct 31, 2018 at 18:27

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.