20190521_eQTL_enrichmentCrossSpecies

workflowr

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.DS_Store
    Ignored:    analysis/20190521_eQTL_CrossSpeciesEnrichment_cache/
    Ignored:    analysis/figure/
    Ignored:    analysis_temp/.DS_Store
    Ignored:    code/.DS_Store
    Ignored:    code/snakemake_workflow/.DS_Store
    Ignored:    data/.DS_Store
    Ignored:    data/PastAnalysesDataToKeep/.DS_Store
    Ignored:    docs/.DS_Store
    Ignored:    docs/assets/.DS_Store

Untracked files:
    Untracked:  analysis/20190617_ExpressionVarianceOfEgenes.Rmd
    Untracked:  docs/figure/20190613_PowerAnalysis.Rmd/

Unstaged changes:
    Modified:   analysis/20190613_PowerAnalysis.Rmd

library(plyr)
library(tidyverse)
library(knitr)
library(data.table)
library(ggpmisc)
library("clusterProfiler")
library("org.Hs.eg.db")
library(gridExtra)
library(enrichplot)

Here the dataset is eQTLs called from a model with the following pre-testing filters/transformations/checks:

• cis-window=250kB
• MAF>10%
• genes tested must have >6 reads in 80% of samples
• lmm with genetic relatedness matrix produced by gemma
• Gene expression is standardized and normalized
• 10PCs added as covariates
• FDR estimated by Storey’s qvalue
• Pvalues well calibrated under a permutated null

First, Read in the data…

eQTLs <- read.table(gzfile("../data/PastAnalysesDataToKeep/20190521_eQTLs_250kB_10MAF.txt.gz"), header=T)
kable(head(eQTLs))

# List of chimp tested genes
ChimpTestedGenes <- rownames(read.table('../output/ExpressionMatrix.un-normalized.txt.gz', header=T, check.names=FALSE, row.names = 1))

ChimpToHumanGeneMap <- read.table("../data/Biomart_export.Hsap.Ptro.orthologs.txt.gz", header=T, sep='\t', stringsAsFactors = F)
kable(head(ChimpToHumanGeneMap))

# Of this ortholog list, how many genes are one2one
table(ChimpToHumanGeneMap$Chimpanzee.homology.type)

ortholog_many2many  ortholog_one2many   ortholog_one2one 
              2278              19917             140351 

OneToOneMap <- ChimpToHumanGeneMap %>%
  filter(Chimpanzee.homology.type=="ortholog_one2one")

# Read gtex heart egene list
# Only consider those that were tested in both species and are one2one orthologs
GtexHeartEgenes <- read.table("../data/Heart_Left_Ventricle.v7.egenes.txt.gz", header=T, sep='\t', stringsAsFactors = F) %>%
  mutate(gene_id_stable = gsub(".\\d+$","",gene_id)) %>%
  filter(gene_id_stable %in% OneToOneMap$Gene.stable.ID) %>%
  mutate(chimp_id = plyr::mapvalues(gene_id_stable, OneToOneMap$Gene.stable.ID,  OneToOneMap$Chimpanzee.gene.stable.ID, warn_missing = F)) %>%
  filter(chimp_id %in% ChimpTestedGenes)

ChimpToHuman.ID <- function(Chimp.ID){
  #function to convert chimp ensembl to human ensembl gene ids
  return(
    plyr::mapvalues(Chimp.ID, OneToOneMap$Chimpanzee.gene.stable.ID, OneToOneMap$Gene.stable.ID, warn_missing = F)
  )}

Now compare with GTEx by making 2x2 contigency table (eGene/not-eGene in Chimp/human). The odds ratio from this table is symetrical.

HumanFDR <- 0.1
ChimpFDR <- 0.1

#Get chimp eQTLs
Chimp_eQTLs <- eQTLs %>%
  filter(qvalue<ChimpFDR)

# Count chimp eGenes
length(unique(Chimp_eQTLs$gene))

[1] 336

# Count human eGenes
length(GtexHeartEgenes %>% filter(qval< HumanFDR) %>% pull(chimp_id))

[1] 5410

# Count number genes tested in both species (already filtered for 1to1 orthologs)
length(GtexHeartEgenes$gene_id_stable)

[1] 11586

The number of human eGenes is huge (about half of all tested genes) and GTEx over-powered compared to chimp. With huge power, everything is an eGene and the eGene classification becomes devoid of meaningful information. So I will play with different ways to classify human eGenes.

#Change FDR thresholds or take top N eGenes by qvalue
HumanTopN <- 600
HumanFDR <- 0.1
ChimpFDR <- 0.1

# Filter human eGenes by qval threshold
HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% pull(chimp_id)

# Filter human eGenes by topN qval
# HumanSigGenes <- GtexHeartEgenes %>% top_n(-HumanTopN, qval) %>% pull(chimp_id)

# Filter human eGeness by qval threshold then topN betas 
# HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% top_n(1000, abs(slope)) %>% pull(chimp_id)

HumanNonSigGenes <- GtexHeartEgenes %>%
  filter(!chimp_id %in% HumanSigGenes) %>%
  pull(chimp_id)

ChimpSigGenes <- GtexHeartEgenes %>%
  filter(chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)

ChimpNonSigGenes <- GtexHeartEgenes %>%
  filter(! chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)



ContigencyTable <- matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),
                              length(intersect(HumanSigGenes,ChimpNonSigGenes)),
                              length(intersect(ChimpSigGenes,HumanNonSigGenes)),
                              length(intersect(ChimpNonSigGenes,HumanNonSigGenes))), 
                           nrow = 2)

rownames(ContigencyTable) <- c("Chimp eGene", "Not Chimp eGene")
colnames(ContigencyTable) <- c("Human eGene", "Not human eGene")

#what is qval threshold for human eGene classification in this contigency table
print(GtexHeartEgenes %>% top_n(-HumanTopN, qval) %>% top_n(1, qval) %>% pull(qval))

[1] 5.83223e-12

#Contigency table of one to one orthologs tested in both chimps and humans of whether significant in humans, or chimps, or both, or neither
ContigencyTable

                Human eGene Not human eGene
Chimp eGene             145             135
Not Chimp eGene        5265            6041

#One-sided Fisher test for greater overlap than expected by chance
fisher.test(ContigencyTable, alternative="greater")

    Fisher's Exact Test for Count Data

data:  ContigencyTable
p-value = 0.04779
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
 1.002643      Inf
sample estimates:
odds ratio 
  1.232354 

The above contingency table and one-sided fisher test indicated a greater-than-chance overlap between the sets of eGenes in chimp and human.

A few plots might give me a better intuition for thinking about the data…

# Plot betas vs qval for GTEx for qval<0.1
GtexHeartEgenes %>%
  filter(qval<0.1) %>%
  ggplot(aes(x=slope, y=-log10(qval))) +
         geom_point(alpha=0.05)

Past versions of random-plots-1.png

HumanAndChimpEffectSizes <- Chimp_eQTLs %>%
    group_by(gene) %>%
    dplyr::slice(which.min(qvalue)) %>% 
    ungroup() %>%
    left_join((GtexHeartEgenes %>% filter(qval<0.1)), by=c("gene"="chimp_id")) %>%
    dplyr::select(gene, slope, beta)

Warning: Column `gene`/`chimp_id` joining factor and character vector,
coercing into character vector

# For shared eGenes, are scatter plot of absolute value of betas
ggplot(HumanAndChimpEffectSizes, aes(x=log10(abs(slope)), y=log10(abs(beta)))) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw() +
  xlim(c(-1,1)) +
  ylim(c(-1,1)) +
  xlab("Human |beta|") +
  ylab("Chimp |beta|")

Warning: Removed 191 rows containing non-finite values (stat_smooth).

Warning: Removed 191 rows containing missing values (geom_point).

Past versions of random-plots-2.png

cor.test(HumanAndChimpEffectSizes$slope, HumanAndChimpEffectSizes$beta,method='spearman')

    Spearman's rank correlation rho

data:  HumanAndChimpEffectSizes$slope and HumanAndChimpEffectSizes$beta
S = 476920, p-value = 0.4633
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
0.06132499 

Effect sizes for cross-species eGenes appear uncorrelated. This might at first appear at odds with the greater-than-chance overlap for the same set of eGenes across sepecies. However, one explanation is that the above effect sizes analysis didn’t consider non eGenes for which the null of effect-size=0 could not be rejected.

OK, so now I will perform the same Odds ratio estimate and Fisher exact test for different thresholds/quantiles of the human data ordered by significance.

ChimpFDR <- 0.1
HumanFDR <- 0.1
Chimp_eQTLs <- eQTLs %>%
  filter(qvalue<ChimpFDR)

ChimpSigGenes <- GtexHeartEgenes %>%
  filter(chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)

ChimpNonSigGenes <- GtexHeartEgenes %>%
  filter(! chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)

LengthOut=200
TopHumanEgeneCountList <- round(seq(100,5410, length.out=LengthOut))
OddsRatiosList <- rep(0,LengthOut)
CI.lower <- rep(0,LengthOut)
CI.upper <- rep(0,LengthOut)

for (i in seq_along(TopHumanEgeneCountList)){
  HumanTopN <- TopHumanEgeneCountList[i]


# Filter human eGenes by qval threshold
# HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% pull(chimp_id)

# Filter human eGenes by topN qval
HumanSigGenes <- GtexHeartEgenes %>% top_n(-HumanTopN, qval) %>% pull(chimp_id)

# Filter human eGeness by qval threshold then topN betas 
# HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% top_n(1000, abs(slope)) %>% pull(chimp_id)

HumanNonSigGenes <- GtexHeartEgenes %>%
  filter(!chimp_id %in% HumanSigGenes) %>%
  pull(chimp_id)

ContigencyTable <- matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),
                              length(intersect(HumanSigGenes,ChimpNonSigGenes)),
                              length(intersect(ChimpSigGenes,HumanNonSigGenes)),
                              length(intersect(ChimpNonSigGenes,HumanNonSigGenes))), 
                           nrow = 2)

#Contigency table of one to one orthologs tested in both chimps and humans of whether significant in humans, or chimps, or both, or neither
ContigencyTable
A<-fisher.test(ContigencyTable, conf.int=T)
OddsRatiosList[i]<- as.numeric(A$estimate)
CI.lower[i] <- as.numeric(A$conf.int[1])
CI.upper[i] <- as.numeric(A$conf.int[2])
}
HeartOddsRatios <- data.frame(TopHumanEgeneCountList=TopHumanEgeneCountList, OddsRatio=OddsRatiosList, GTEx.Tissue="Heart_LeftVentricle", CI.lower=CI.lower, CI.upper=CI.upper)

ggplot(HeartOddsRatios, aes(x=TopHumanEgeneCountList, y=OddsRatiosList)) +
  geom_point() +
  geom_ribbon(aes(ymin=CI.lower, ymax=CI.upper), alpha=0.2) +
  theme_bw() +
  ylab("Odds Ratio") +
  xlab("Considering only top X human eGenes") +
  geom_hline(yintercept=1, linetype="dashed")

Past versions of CompareWithGtexEgenes_byquantiles-1.png

Now do the same for a couple of GTEx tissues

GTEx.eGene.Files <- c("../data/Heart_Left_Ventricle.v7.egenes.txt.gz",
                      "../data/Adipose_Subcutaneous.v7.egenes.txt.gz",
                      # "../data/Cells_EBV-transformed_lymphocytes.v7.egenes.txt.gz",
                      "../data/Thyroid.v7.egenes.txt.gz",
                      "../data/Lung.v7.egenes.txt.gz")

GTEx.tissue.labels <- c("Heart_Left_Ventricle",
                        "Adipose_Subcutaneous",
                        # "LCL",
                        "Thyroid",
                        "Lung")

OddsRatios.rbinded <- data.frame()
                      
for (j in seq_along(GTEx.eGene.Files)){
  paste0("processing")
  GtexEgenes <- read.table(GTEx.eGene.Files[j], header=T, sep='\t', stringsAsFactors = F) %>%
  mutate(gene_id_stable = gsub(".\\d+$","",gene_id)) %>%
  filter(gene_id_stable %in% OneToOneMap$Gene.stable.ID) %>%
  mutate(chimp_id = plyr::mapvalues(gene_id_stable, OneToOneMap$Gene.stable.ID,  OneToOneMap$Chimpanzee.gene.stable.ID, warn_missing = F)) %>%
  filter(chimp_id %in% ChimpTestedGenes)

LengthOut=200
TopHumanEgeneCountList <- round(seq(100,5410, length.out=LengthOut))
OddsRatiosList <- rep(0,LengthOut)
CI.lower <- rep(0,LengthOut)
CI.upper <- rep(0,LengthOut)

for (i in seq_along(TopHumanEgeneCountList)){
  HumanTopN <- TopHumanEgeneCountList[i]


  # Filter human eGenes by qval threshold
  # HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% pull(chimp_id)
  
  # Filter human eGenes by topN qval
  HumanSigGenes <- GtexEgenes %>% top_n(-HumanTopN, qval) %>% pull(chimp_id)
  
  HumanNonSigGenes <- GtexEgenes %>%
    filter(!chimp_id %in% HumanSigGenes) %>%
    pull(chimp_id)
  
  ContigencyTable <- matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),
                                length(intersect(HumanSigGenes,ChimpNonSigGenes)),
                                length(intersect(ChimpSigGenes,HumanNonSigGenes)),
                                length(intersect(ChimpNonSigGenes,HumanNonSigGenes))), 
                             nrow = 2)
  
  #Contigency table of one to one orthologs tested in both chimps and humans of whether significant in humans, or chimps, or both, or neither
  A<-fisher.test(ContigencyTable, conf.int=T)
  OddsRatiosList[i]<- as.numeric(A$estimate)
  CI.lower[i] <- as.numeric(A$conf.int[1])
  CI.upper[i] <- as.numeric(A$conf.int[2])
}

OddsRatios <- data.frame(HumanEgeneThreshold=TopHumanEgeneCountList, OR=OddsRatiosList, GTEx.Tissue=GTEx.tissue.labels[j], CI.lower=CI.lower, CI.upper=CI.upper)

OddsRatios.rbinded <- rbind(OddsRatios.rbinded, OddsRatios)
  
}

ggplot(OddsRatios.rbinded, aes(x=HumanEgeneThreshold, y=OR, color=GTEx.Tissue)) +
  geom_line(size=2) +
  geom_ribbon(aes(ymin=CI.lower, ymax=CI.upper, fill = GTEx.Tissue), alpha=0.1,linetype = 3) +
  theme_bw() +
  ylab("Odds Ratio") +
  xlab("Considering only top X human eGenes") +
  geom_hline(yintercept=1, linetype="dashed") +
  ylim(c(0,4)) +
  theme(legend.position="top")

Past versions of CompareWithGtexEgenes_OtherTissue-1.png

Try one alternative way to present these odds ratios as a function of eGene stringency:

# Get qval ordered ordered list of Chimp genes (that were tested in both)
Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    arrange(qvalue) %>%
    as.data.frame()

# Get qval ordered list of human genes 
Human_OrderedGenes <- GtexHeartEgenes %>%
  arrange(qval)


LengthOut=60

ChimpEgeneCount <- Chimp_OrderedGenes %>% filter(qvalue<0.2) %>% pull(qvalue) %>% length()
HumanEgeneCount <- Human_OrderedGenes %>% filter(qval<0.2) %>% pull(qval) %>% length()
TopHumanEgeneCountList <- round(seq(1,HumanEgeneCount, by=50))[-1]
TopChimpEgeneCountList <- round(seq(1,ChimpEgeneCount, by=25))[-1]

OddsRatioMatrix <- matrix(nrow = length(TopHumanEgeneCountList), ncol = length(TopChimpEgeneCountList))

for (i in seq_along(TopHumanEgeneCountList)){
  HumanTopN <- TopHumanEgeneCountList[i]

# Filter human eGenes by topN qval
HumanSigGenes <- Human_OrderedGenes %>% top_n(-HumanTopN, qval) %>% pull(chimp_id)

HumanNonSigGenes <- Human_OrderedGenes %>%
  filter(!chimp_id %in% HumanSigGenes) %>%
  pull(chimp_id)

for (j in seq_along(TopChimpEgeneCountList)){
  TopChimpN <- TopChimpEgeneCountList[j]
  
  ChimpSigGenes <- Chimp_OrderedGenes %>% top_n(-TopChimpN, qvalue) %>% pull(gene)
  ChimpNonSigGenes <- GtexHeartEgenes %>%
    filter(!chimp_id %in% ChimpSigGenes) %>%
    pull(chimp_id)
  
  ContigencyTable <- matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),
                              length(intersect(HumanSigGenes,ChimpNonSigGenes)),
                              length(intersect(ChimpSigGenes,HumanNonSigGenes)),
                              length(intersect(ChimpNonSigGenes,HumanNonSigGenes))), 
                           nrow = 2)
  A<-fisher.test(ContigencyTable, conf.int=F)
  OddsRatioMatrix[i,j] <- (as.numeric(A$estimate))
}
}

row.names(OddsRatioMatrix) <- TopHumanEgeneCountList
colnames(OddsRatioMatrix) <- TopChimpEgeneCountList

melt(OddsRatioMatrix) %>%
  dplyr::rename(OddsRatio=value) %>%
ggplot(aes(x=as.numeric(Var1),y=as.numeric(Var2), fill=OddsRatio)) +
  geom_raster(interpolate=F) +
  xlab("Top N human eGenes") +
  ylab("Top N chimp eGenes") +
  scale_y_continuous(expand = c(0,0)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_fill_distiller(palette = "Spectral") +
  # coord_fixed(ratio = 1, xlim = NULL, ylim = NULL, expand = F, clip = "on") +
  
  geom_vline(xintercept=7555, colour="grey") +
  geom_text(aes(x=7555, label="FDR=0.2", y=400), angle=90, vjust = 0, color="grey", size=4)+
  geom_vline(xintercept=5410, colour="grey") +
  geom_text(aes(x=5410, label="FDR=0.1", y=400), angle=90, vjust = 0, color="grey", size=4)+
  geom_vline(xintercept=4240, colour="grey") +
  geom_text(aes(x=4240, label="FDR=0.05", y=400), angle=90, vjust = 0, color="grey", size=4)+
  
    geom_hline(yintercept=731, colour="grey") +
  geom_text(aes(x=2000, label="FDR=0.2", y=731-50), vjust = 0, color="grey", size=4)+
  geom_hline(yintercept=280, colour="grey") +
  geom_text(aes(x=2000, label="FDR=0.1", y=280-50), vjust = 0, color="grey", size=4)+
  geom_hline(yintercept=98, colour="grey") +
  geom_text(aes(x=2000, label="FDR=0.05", y=98-50), vjust = 0, color="grey", size=4)+
  theme_bw() +
  theme(legend.position="bottom")

Past versions of unnamed-chunk-1-1.png

Gene ontology overlap analysis of genes that are Chimp eGenes but not human eGenes. The interpretation is that this is a somewhat prioritized list of candidate genes under stabilizing selection in the human lineage.

#Change FDR thresholds or take top N eGenes by qvalue
HumanFDR <- 0.1
ChimpFDR <- 0.1
HumanTopN <- 600

Chimp_eQTLs <- eQTLs %>%
  filter(qvalue<ChimpFDR)

# Filter human eGenes by qval threshold
HumanSigGenes <- GtexHeartEgenes %>% filter(qval<HumanFDR) %>% pull(chimp_id)
# HumanSigGenes <- GtexHeartEgenes %>% top_n(-HumanTopN, qval) %>% pull(chimp_id)


HumanNonSigGenes <- GtexHeartEgenes %>%
  filter(!chimp_id %in% HumanSigGenes) %>%
  pull(chimp_id)

ChimpSigGenes <- GtexHeartEgenes %>%
  filter(chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)

ChimpNonSigGenes <- GtexHeartEgenes %>%
  filter(! chimp_id %in% Chimp_eQTLs$gene) %>%
  pull(chimp_id)



ContigencyTable <- matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),
                              length(intersect(HumanSigGenes,ChimpNonSigGenes)),
                              length(intersect(ChimpSigGenes,HumanNonSigGenes)),
                              length(intersect(ChimpNonSigGenes,HumanNonSigGenes))), 
                           nrow = 2)

rownames(ContigencyTable) <- c("Chimp eGene", "Not Chimp eGene")
colnames(ContigencyTable) <- c("Human eGene", "Not human eGene")


#Contigency table of one to one orthologs tested in both chimps and humans of whether significant in humans, or chimps, or both, or neither
ContigencyTable

                Human eGene Not human eGene
Chimp eGene             145             135
Not Chimp eGene        5265            6041

Now that we have a contingency table of eGenes between chimp and human, lets do gene set enrichment analysis. This analysis uses clusterProfiler package to look for gene set enrichment in a foreground set of genes over a background set. For each foreground/background comparison, there will be a gene ontology enrichment for the GO subcategories of - biological process (bp) - cellular component (cc) - molecular function (mf)

Perhaps first most intuitive comparison to make is:

Foreground = Chimp eGene & not human eGene

background = All genes tested in both species

foreground = ChimpToHuman.ID(intersect(ChimpSigGenes,HumanNonSigGenes))
background = ChimpToHuman.ID(GtexHeartEgenes$chimp_id)
length(foreground)

[1] 135

length(background)

[1] 11586

ego.mf <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "MF", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.mf)))

ego.bp <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "BP", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.bp)))

ego.cc <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "CC", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.cc)))

Nothing enriched by this comparison of gene sets. Perhaps a more appropriate foreground/background set is:

foreground = chimp specific eGenes

background = all chimp eGenes

foreground = ChimpToHuman.ID(intersect(ChimpSigGenes,HumanNonSigGenes))
background = ChimpToHuman.ID(ChimpSigGenes)
length(foreground)

[1] 135

length(background)

[1] 280

ego.mf <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "MF", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.mf)))

ego.bp <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "BP", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.bp)))

ego.cc <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "CC", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.cc)))

dotplot(ego.mf, font.size=8)

Past versions of “GO-2, -1.png”

dotplot(ego.bp, font.size=8)

Past versions of “GO-2, -2.png”

dotplot(ego.cc, font.size=8)

Past versions of “GO-2, -3.png”

Compare those that are shared species eGenes, with the interpretation that they would be enriched for things not under stabilizing selection in either, perhaps things like relatively useless pseudogenes.

foreground = shared eGenes

background = all genes tested in both species

foreground = ChimpToHuman.ID(intersect(ChimpSigGenes,HumanSigGenes))
background = ChimpToHuman.ID(GtexHeartEgenes$chimp_id)
length(foreground)

[1] 145

length(background)

[1] 11586

# foreground = ChimpToHuman.ID(intersect(ChimpSigGenes,(GtexHeartEgenes %>% top_n(-600, qval) %>% pull(chimp_id))))
# Could also use the more stringent intersection of chimp eGenes and top600 human eGenes. You get mostly the sae GO terms, including MHC complex

ego.mf <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "MF", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.mf)))

ego.bp <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "BP", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.bp)))

ego.cc <- enrichGO(gene = foreground,
                 universe = background,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont = "CC", 
                 pAdjustMethod = "BH", 
                 pvalueCutoff  = 0.01,
                 qvalueCutoff  = 0.05)
kable(head(as.data.frame(ego.cc)))

dotplot(ego.mf, font.size=8)

Past versions of GO-3-1.png

dotplot(ego.bp, font.size=8)

Past versions of GO-3-2.png

dotplot(ego.cc, font.size=8)

Past versions of GO-3-3.png

Interestingly, the MHC complex comes out of this analysis as enriched among the genes that are eGenes in both species, consistent with these genes being rapidly evolving targets of selection in between human/chimp, with high diversity in both species.

A complementart analysis to GO overlap is gene set enrichment analysis. Here I will perform GSEA on the chimp qvalue orded list of chimp/human shared eGenes and chimp specific eGenes based on FDR<0.1 in the human dataset…

Actually, maybe a more sensible thing to do is consider the human qvalue ordered list and then separate based on if chimp eGene or not. This is more similar to what the previous odds-ratio plots are.

HumanOrderedGenes <- GtexHeartEgenes %>%
  arrange(qval) %>%
  filter(chimp_id %in% ChimpSigGenes)
GeneList <- log10(HumanOrderedGenes$qval) * -1
names(GeneList) <- as.character(ChimpToHuman.ID(HumanOrderedGenes$chimp_id))

gsego.cc.shared <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 # pvalueCutoff=0.8,
                 ont = "CC")
kable(head(as.data.frame(gsego.cc.shared)))

dotplot(gsego.cc.shared, font.size=8, showCategory=15)
gseaplot(gsego.cc.shared, geneSetID = "GO:0042611", pvalue_table = TRUE)

That was still a work in progress that I may come back to one day… For now I will perform GSEA on ranked gene lists based on chimp FDR.

First, create chimp-qval-ordered gene list for genes that are human eGenes I expect this to return similar results as gene overlap analysis with forground as shared eGenes and background as all all genes tested. So this should return MHC complex among other things.

Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    arrange(qvalue) %>%
    filter(gene %in% HumanSigGenes)
GeneList <- log10(Chimp_OrderedGenes$qvalue) * -1
names(GeneList) <- as.character(ChimpToHuman.ID(Chimp_OrderedGenes$gene))



gsego.cc.shared <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 # pvalueCutoff=0.8,
                 ont = "CC",
                 nPerm = 10000)

preparing geneSet collections...

GSEA analysis...

Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, : There are ties in the preranked stats (38.71% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

leading edge analysis...

done...

kable(head(as.data.frame(gsego.cc.shared)))

dotplot(gsego.cc.shared, font.size=8, showCategory=15)

Past versions of unnamed-chunk-3-1.png

gseaplot2(gsego.cc.shared, geneSetID = c("GO:0042611", "GO:0042613", "GO:0032588"), title="Shared eGenes ordered by chimp -log(qvalue)", pvalue_table = TRUE)

Past versions of unnamed-chunk-3-2.png

Ok that does return MHC as expect. Now I will create chimp-qval-ordered gene list for genes that are not human eGenes. I expect this to not return MHC. Instead,I expect this to return similar results as gene overlap analysis with forground as chimp-specific eGenes and background as all chimp-eGenes (nucleic acid binding, etc).

Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    filter(qvalue<0.2) %>%
    arrange(qvalue) %>%
    filter(!gene %in% HumanSigGenes)
GeneList <- log10(Chimp_OrderedGenes$qvalue) * -1
names(GeneList) <- as.character(ChimpToHuman.ID(Chimp_OrderedGenes$gene))

gsego.cc.chimpspecific <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 pvalueCutoff=0.8,
                 ont = "BP")

preparing geneSet collections...

GSEA analysis...

Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, : There are ties in the preranked stats (28.95% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

no term enriched under specific pvalueCutoff...

kable(head(as.data.frame(gsego.cc.chimpspecific)))

Ok that returned nothing, even at more relaxed pvalue thresholds (not shown) I can see the top terms do not have anything to do with nucleic acid binding like the previous gene overlap enrichment analysis results that compared chimp-specific-eGene (foreground) to all chimp eGenes (universe background).

By partitioning the GSEA analyses into “human eGene” ranked gene list and “not human eGene” ranked gene list, what I was really trying to get at is the semi-quantitative difference in eGene character. Essentially these are paired analyses, though GSEA does not have any natural way to run a paired analysis on two gene sets. Perhaps a better way to do this is to make a single gene list based on the difference in eGene ranks between the two species. Since GSEA is based on non-parametric rank orders, I think this is still somewhat reasonable to do. For example, if qvalues are ranked lowest to highest in both species (rank=1 being most significant), and I subtract the chimp ranking from the human ranking, I am left with the more chimp specific eGenes with a higher number. I will implement this idea below:

Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    left_join(GtexHeartEgenes, by=c("gene"="chimp_id")) %>%
    dplyr::select(gene, qvalue, qval) %>% as.data.frame() %>%
    mutate(ChimpRank = dense_rank(qvalue)) %>%
    mutate(HumanRank = dense_rank(qval)) %>%
    mutate(RankDifference = HumanRank-ChimpRank) %>%
    arrange(desc(RankDifference))

Warning: Column `gene`/`chimp_id` joining factor and character vector,
coercing into character vector

GeneList <- Chimp_OrderedGenes$RankDifference
names(GeneList) <- as.character(ChimpToHuman.ID(Chimp_OrderedGenes$gene))
    
gsego.cc.chimpspecific <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 pvalueCutoff=0.8,
                 ont = "BP")

preparing geneSet collections...

GSEA analysis...

Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, : There are ties in the preranked stats (35.5% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

leading edge analysis...

done...

kable(head(as.data.frame(gsego.cc.chimpspecific)))

# gseaplot(gsego.cc.chimpspecific, geneSetID = "GO:0001501", pvalue_table = TRUE)
# dotplot(gsego.cc.chimpspecific, font.size=8, showCategory=15)

Nothing is significant by any notable treshold by this method…

Perhaps I need to limit the gene set to genes which were called as eGenes in chimps, since many too many chimp non-eGenes may contribute to noise.

Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    left_join(GtexHeartEgenes, by=c("gene"="chimp_id")) %>%
    dplyr::select(gene, qvalue, qval) %>% as.data.frame() %>%
    mutate(ChimpRank = dense_rank(qvalue)) %>%
    mutate(HumanRank = dense_rank(qval)) %>%
    mutate(RankDifference = HumanRank-ChimpRank) %>%
    filter(qvalue <0.1) %>%
    arrange(desc(RankDifference))

Warning: Column `gene`/`chimp_id` joining factor and character vector,
coercing into character vector

GeneList <- Chimp_OrderedGenes$RankDifference
names(GeneList) <- as.character(ChimpToHuman.ID(Chimp_OrderedGenes$gene))
    
gsego.cc.chimpspecific <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 ont = "BP",
                 nPerm = 100000)

preparing geneSet collections...

GSEA analysis...

Warning in fgsea(pathways = geneSets, stats = geneList, nperm = nPerm, minSize = minGSSize, : There are ties in the preranked stats (2.14% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

leading edge analysis...

done...

kable(head(as.data.frame(gsego.cc.chimpspecific), 2))

gseaplot2(gsego.cc.chimpspecific, geneSetID = c("GO:0006366","GO:0044260"), title="Chimp eGenes ordered by difference in species rank", pvalue_table = TRUE)

Past versions of unnamed-chunk-6-1.png

dotplot(gsego.cc.chimpspecific, font.size=8, showCategory=38)

Past versions of unnamed-chunk-6-2.png

Alright. Great. that analysis recapitulated the gene overlap analysis of [foreground=Chimp-specific-eGenes; background=all-chimp-eGenes] in the sense that terms like ‘gene transcription’ pop up.

Let’s look at some of the genes in this GO category that are at the top of this ranked metric (have more eGene character in chimp than human)

gsea.results <- as.data.frame(gsego.cc.chimpspecific)
PolII.genes.in.set <- gsea.results %>% 
  filter(Description=="transcription by RNA polymerase II") %>% 
  pull(core_enrichment) %>%
  strsplit("/") %>% unlist() %>%
  mapIds(org.Hs.eg.db, keys=., column="SYMBOL", keytype="ENSEMBL", multiVals="first")

'select()' returned 1:1 mapping between keys and columns

PolII.genes.in.set

ENSG00000173575 ENSG00000102974 ENSG00000164758 ENSG00000196150 
         "CHD2"          "CTCF"         "MED30"        "ZNF250" 
ENSG00000152454 ENSG00000094916 ENSG00000078401 ENSG00000157404 
       "ZNF256"          "CBX5"          "EDN1"           "KIT" 
ENSG00000181007 ENSG00000168743 ENSG00000151623 ENSG00000150593 
        "ZFP82"          "NPNT"         "NR3C2"         "PDCD4" 
ENSG00000106261 ENSG00000163602 ENSG00000170325 ENSG00000153187 
      "ZKSCAN1"          "RYBP"        "PRDM10"        "HNRNPU" 
ENSG00000159388 ENSG00000123609 ENSG00000151718 ENSG00000102984 
         "BTG2"           "NMI"          "WWC2"        "ZNF821" 
ENSG00000114503 ENSG00000185122 ENSG00000164342 ENSG00000161057 
        "NCBP2"          "HSF1"          "TLR3"         "PSMC2" 
ENSG00000141646 ENSG00000152409 ENSG00000204519 
        "SMAD4"           "JMY"        "ZNF551" 

CTCF (ENSG00000102974) and HSF1 (ENSG00000185122) sticks out as something interesting and recognizeable genes. Also, quite a few zinc-finger proteins (ZNF821, ZNF551, ZNF250, ZFP82, ZNF256, PRDM10) Let’s see exactly where they are on this ranked list…

# All the top genes by this difference in ranks ordered gene list
GeneSymbolAdded <- Chimp_OrderedGenes %>% 
  mutate(geneSymbol = mapIds(org.Hs.eg.db, keys=ChimpToHuman.ID(gene), column="SYMBOL", keytype="ENSEMBL", multiVals="first") )

'select()' returned 1:many mapping between keys and columns

GeneSymbolAdded %>% head(40) %>% kable()

# The "interesting genes" that are in the PolII regulation of transcription gene set
GeneSymbolAdded %>% 
  filter(geneSymbol %in% PolII.genes.in.set) %>% kable()

CTCF and other zinc fingers are among the top in the list (some top chimp eGenes but not for human). This surprises me considering I think of CTCF as so conserved across large evolutionary distance. I will have to look back at more raw data to reassure myself that these chimp eQTLs are real.

Now let’s try a reciprocal analysis; where we look for what is enriched in human-specific eGenes, (again where the metric to order the gene list is a difference in ranks.)… Still a work-in-progress

Chimp_OrderedGenes <- eQTLs %>%
    group_by(gene) %>% 
    dplyr::slice(which.min(qvalue)) %>%
    filter(gene %in% GtexHeartEgenes$chimp_id) %>%
    filter(qvalue <0.1) %>%
    left_join(GtexHeartEgenes, by=c("gene"="chimp_id")) %>%
    dplyr::select(gene, qvalue, qval) %>% as.data.frame() %>%
    mutate(ChimpRank = dense_rank(qvalue)) %>%
    mutate(HumanRank = dense_rank(qval)) %>%
    mutate(RankDifference = HumanRank-ChimpRank) %>%
    arrange(desc(RankDifference))
GeneList <- Chimp_OrderedGenes$RankDifference
names(GeneList) <- as.character(ChimpToHuman.ID(Chimp_OrderedGenes$gene))
    
gsego.cc.chimpspecific <- gseGO(gene = GeneList,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 maxGSSize = 500,
                 ont = "BP",
                 nPerm = 10000)
kable(head(as.data.frame(gsego.cc.chimpspecific), 15))
gseaplot2(gsego.cc.chimpspecific, geneSetID = c("GO:0006357", "GO:0006355"), title="Chimp eGenes ordered by difference in species rank", pvalue_table = TRUE)

dotplot(gsego.cc.chimpspecific, font.size=8, showCategory=60)

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] enrichplot_1.2.0       gridExtra_2.3          org.Hs.eg.db_3.7.0    
 [4] AnnotationDbi_1.44.0   IRanges_2.16.0         S4Vectors_0.20.1      
 [7] Biobase_2.42.0         BiocGenerics_0.28.0    clusterProfiler_3.10.1
[10] ggpmisc_0.3.1          data.table_1.12.2      knitr_1.23            
[13] forcats_0.4.0          stringr_1.4.0          dplyr_0.8.1           
[16] purrr_0.3.2            readr_1.3.1            tidyr_0.8.3           
[19] tibble_2.1.3           ggplot2_3.1.1          tidyverse_1.2.1       
[22] plyr_1.8.4            

loaded via a namespace (and not attached):
 [1] nlme_3.1-140        fs_1.3.1            lubridate_1.7.4    
 [4] bit64_0.9-7         progress_1.2.2      RColorBrewer_1.1-2 
 [7] httr_1.4.0          UpSetR_1.4.0        rprojroot_1.3-2    
[10] tools_3.5.1         backports_1.1.4     R6_2.4.0           
[13] DBI_1.0.0           lazyeval_0.2.2      colorspace_1.4-1   
[16] withr_2.1.2         prettyunits_1.0.2   tidyselect_0.2.5   
[19] bit_1.1-14          compiler_3.5.1      git2r_0.25.2       
[22] cli_1.1.0           rvest_0.3.4         xml2_1.2.0         
[25] labeling_0.3        triebeard_0.3.0     scales_1.0.0       
[28] ggridges_0.5.1      digest_0.6.19       rmarkdown_1.13     
[31] DOSE_3.8.2          pkgconfig_2.0.2     htmltools_0.3.6    
[34] highr_0.8           rlang_0.3.4         readxl_1.3.1       
[37] rstudioapi_0.10     RSQLite_2.1.1       gridGraphics_0.4-1 
[40] generics_0.0.2      farver_1.1.0        jsonlite_1.6       
[43] BiocParallel_1.16.6 GOSemSim_2.8.0      magrittr_1.5       
[46] ggplotify_0.0.3     GO.db_3.7.0         Matrix_1.2-17      
[49] Rcpp_1.0.1          munsell_0.5.0       viridis_0.5.1      
[52] stringi_1.4.3       whisker_0.3-2       yaml_2.2.0         
[55] ggraph_1.0.2        MASS_7.3-51.4       qvalue_2.14.1      
[58] grid_3.5.1          blob_1.1.1          ggrepel_0.8.1      
[61] DO.db_2.9           crayon_1.3.4        lattice_0.20-38    
[64] cowplot_0.9.4       haven_2.1.0         splines_3.5.1      
[67] hms_0.4.2           pillar_1.4.1        fgsea_1.8.0        
[70] igraph_1.2.4.1      reshape2_1.4.3      fastmatch_1.1-0    
[73] glue_1.3.1          evaluate_0.14       modelr_0.1.4       
[76] urltools_1.7.3      tweenr_1.0.1        cellranger_1.1.0   
[79] gtable_0.3.0        polyclip_1.10-0     assertthat_0.2.1   
[82] xfun_0.7            ggforce_0.2.2       europepmc_0.3      
[85] broom_0.5.2         viridisLite_0.3.0   rvcheck_0.1.3      
[88] memoise_1.1.0       workflowr_1.4.0