How to export in tsv file, when I click on download option it opens a new window with tex.

There is no option for saving as.

Hi Kevin,

I am glad I found your thread, as it is basically what I am trying to do/present my data, but being new to R so I am finding it tricky.

I have Normalised gene expression data (mRNA) and I have managed to use DAVID to cluster the Genes, and have the Enrichment file.

I want to display the mean normalised values rather than significance. i.e A heatmap to display upregulated genes and another heatmap to display downregulated genes in the pathways.

I don't know how to statistically analyse the data using DESeq2, I have calculated the P-Value on Excel, so I only have the Mean of my normalised data set and the P-Value (no base mean, no log 2 fold change, no adjusted P-value), is it possible to generate a heatmap with just the data sets I have? I just want to have simple heatmaps, as I want to see which pathway/s are mostly dysregulated in response to my treatment, and then start to create other graphical representations of interesting pathways/genes.

I have tried to follow your example above, but R gives me an error for rldMatrix (saying it could not find the function). I have managed to write the codes for the rest, and it seems to be fine. Just finding the initial bit (data organisation) tricky.

Thank you for your help. Sehar :)

@Kevin you have a topMatrix object here , but i did clustering and taking some of the clustered gene set and giving that as input to david so what is my top matrix as im doing wgcna so there is no differential expression analysis involved so the only thing i have gene list from modules ..any suggestion

Very useful and a great way to look that data will multi feature. Will give it a try. Great resource to be bookmarked.

Hi@Kevin now since i have to do pairwise comparison im back it im doing every thing as you have mentioned after running the for loop i get this
when i see annGSEA i get this "disulfide bond Disulfide bond glycosylation site:N-linked (GlcNAc...) Glycoprotein signal peptide GO:0005615~extracellular space Polymorphism Signal hsa05310:Asthma topological domain:Extracellular GO:0005886~plasma membrane Receptor GO:0005887~integral component of plasma membrane Membrane as my colnames and my first column being my gene names so all my ontology columns are zero im not sure what im doing wrong..

log2FCcutoff <- 3
    BHcutoff <- 0.05
    sigGeneList <- subset(topTable, abs(log2FoldChange)>=log2FCcutoff & padj<=BHcutoff)[,1]




head(sigGeneList)


DAVIDfile <- "MyEnrichment.tsv"
DAVID <- read.table(DAVIDfile, sep="\t", header=TRUE)
colnames(DAVID)


enrichBcutoff <- 0.05
DAVID <- subset(DAVID, Benjamini<enrichBcutoff)
DAVID <- DAVID[,c(1,2,3,6)]



annGSEA <- data.frame(row.names=sigGeneList)
for (j in 1:length(sigGeneList))
{
  #Create a matching pattern to ensure genes match exactly
  #   ^GENE,  --> Match at beginning of matching string
  #   , GENE$ --> Match at end of matching string
  #    GENE,  --> Match between first and last gene in matching string
  gene <- sigGeneList[j]
  #pattern <- paste("^", gene, ", |, ", gene, "$| ", gene, ",", sep="")
  pattern <- paste("^", gene, ", |, ", gene, "$| ", gene, ",", sep="")

  for (k in 1:nrow(DAVID))
  {
    if (any(grep(pattern, DAVID$Genes[k])))
    {
      annGSEA[j,k] <- 1
    }
    else
    {
      annGSEA[j,k] <- 0
    }
  }
}
colnames(annGSEA) <- DAVID[,2]

#Remove terms with no overlapping genes
annGSEA <- annGSEA[,apply(annGSEA, 2, mean)!=0]

#Remove genes with no overlapping terms
annGSEA <- annGSEA[apply(annGSEA, 1, mean)!=0,]

this the code exact code im running ,my "MyEnrichment.tsv" contains the DAVID output which i get after giving gene list for enrichment with cutoff range .so im not sure what im doing wrong

Hello @kevin, Thank you for the code ! It is really helpful My question is every time I plot new dataset with this code, colors for the DEG sig score are keep changing. How can I fix the color of DEG significance color bar (as well as GO term score bar), Thanks in advance,

Hello Kevin Blighe, I have one question. I did Functional Annotation Clustering and i got 100 clusters. But on the annotation results Ensemble ID is coming . Is there any way to replace those ensemble ID to gene name from functional annotation clustering output

Many thanks Kevin for this great tutorial. I went through all the step (just for up-regulated genes) and the image that I got at the end is like this. Why at the bottom of the right side there are some muddled text. How can I remove them? I tried but I couldn't. I really appreciate your time and help.

Hi I am quite new with coding and cant seem to be able to solve/understand this error

Error in grepl(pattern, DAVID$Genes[k]) : invalid regular expression '^c(220.7243738, 65.54562172,

I've copied and pasted the code you have provided so I've no idea what went wrong.

annGSEA <- data.frame(row.names=sigGeneList)
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
    }
  }
}

It would be great if you could share the updated version of the tutorial for all since I'm sure there will be many beginners like me. I'm assuming the output you are referring to is before enrichment cutoff.

> str(DAVID)
'data.frame':   495 obs. of  13 variables:
 $ Category       : Factor w/ 13 levels "BBID","BIOCARTA",..: 12 4 12 13 4 12 12 8 6 13 ...
 $ Term           : Factor w/ 495 levels "109.Chemokine_families",..: 30 94 33 469 95 51 462 332 78 32 ...
 $ Count          : int  35 42 120 116 42 139 84 36 28 104 ...
 $ %              : num  10.8 13 37.2 35.9 13 ...
 $ PValue         : num  4.30e-27 5.77e-22 2.22e-21 2.67e-20 2.93e-20 ...
 $ Genes          : Factor w/ 392 levels "ADAMTS14, RSPO3, UNC5C, THBS1",..: 14 16 270 268 71 314 235 15 373 269 ...
 $ List Total     : int  304 279 304 297 279 304 304 159 276 297 ...
 $ Pop Hits       : int  190 379 3434 3346 421 4551 1965 243 176 2917 ...
 $ Pop Total      : int  20581 16792 20581 20063 16792 20581 20581 6879 16881 20063 ...
 $ Fold Enrichment: num  12.47 6.67 2.37 2.34 6 ...
 $ Bonferroni     : num  1.20e-24 9.51e-19 6.21e-19 2.57e-17 4.84e-17 ...
 $ Benjamini      : num  1.20e-24 9.51e-19 3.11e-19 2.57e-17 2.42e-17 ...
 $ FDR            : num  5.69e-24 9.70e-19 2.93e-18 4.20e-17 4.93e-17 ...

> DAVID[1:5,1:5]
          Category                             Term Count        %   PValue
1      UP_KEYWORDS                         Cytokine    35 10.83591 4.30e-27
2 GOTERM_BP_DIRECT GO:0006954~inflammatory response    42 13.00310 5.77e-22
3      UP_KEYWORDS                   Disulfide bond   120 37.15170 2.22e-21
4   UP_SEQ_FEATURE                   signal peptide   116 35.91331 2.67e-20
5 GOTERM_BP_DIRECT       GO:0006955~immune response    42 13.00310 2.93e-20

Hi Kevin, thanks for this great tutorial. However, can you please elaborate on how did you calculated the gene score? Thanks

Hello @Kevin how can i download DAVID enrichment file in .david format ?

can we do the same with multiple samples instead of genes? Like samples as column and pathway as row?

This is a really great piece of code, Kevin, thank you! I had to modify the first few lines from the files on github to get things going, so sharing here in case the holdup of unique row names is catching anyone else new to R. Reminder to everyone that to pull a csv from github you just need to open the raw file, save as excel and excel does a pretty good job with the rest, don't let troubleshooting some small problems with the csv file hold you back from experimenting with this type of heatmap.

##Organize the results table

topTable <- read.table('Post299161_DEresults.csv', header = TRUE, stringsAsFactors = FALSE, sep = ',')
row.names(topTable) <- paste(make.unique(topTable[,1]))
topTable <- topTable[,-1]
head(topTable)

The bar for normalised enrichment score is overlaying the enriched terms. what could be the reason for this error?