Plot heatmap in many styles

pi1

First let’s read in data for Pi1 method

FilesChunks <- paste0("../code/scratch/PairwisePi1Traits.P.", 1:10, ".txt.gz")


dat <- lapply(FilesChunks, fread, sep='\t') %>%
  bind_rows()

RecodeDat <- read_tsv("../data/Phenotypes_recode_for_Plotting.txt")

RecodeVec <- RecodeDat %>%
  select(PC, ShorterAlias) %>%
  deframe()

RecodeIncludePCs <- RecodeDat %>%
  filter(Include) %>%
  pull(PC)

colnames(dat)

 [1] "PC1"                      "P1"                      
 [3] "GeneLocus"                "p_permutation.x"         
 [5] "singletrait_topvar.x"     "singletrait_topvar_chr.x"
 [7] "singletrait_topvar_pos.x" "FDR.x"                   
 [9] "PC2"                      "P2"                      
[11] "p_permutation.y"          "singletrait_topvar.y"    
[13] "singletrait_topvar_chr.y" "singletrait_topvar_pos.y"
[15] "FDR.y"                    "trait.x.p.in.y"          

I also need to simulate the null distribution to deal with the one to many problem…

MaxSampleSizeToCreateANull <- 100
NSamplesToEstimateDistribution <- 10000
NullSimulatedTestStats <- matrix(nrow=MaxSampleSizeToCreateANull, ncol=NSamplesToEstimateDistribution)
rownames(NullSimulatedTestStats) <- paste0("runif_samplesize", 1:MaxSampleSizeToCreateANull)
colnames(NullSimulatedTestStats) <- paste0("Sample_", 1:NSamplesToEstimateDistribution)


for (i in 1:MaxSampleSizeToCreateANull){
  SampleSizeFromUniform <- i
  SampledDat <- matrix(runif(SampleSizeFromUniform*NSamplesToEstimateDistribution), nrow=NSamplesToEstimateDistribution)
  SampleDatNullTestStatistics <- -log10(apply(SampledDat, 1, min))
  NullSimulatedTestStats[i,] <- SampleDatNullTestStatistics
}


NullSimulatedTestStats %>%
  as.data.frame() %>%
  rownames_to_column("runif_samplesize") %>%
  slice(1:20) %>%
  # mutate(runif_samplesize = as.numeric(str_remove(runif_samplesize, "runif_samplesize"))) %>%
  gather(key="Sample", value="value", -runif_samplesize) %>%
  ggplot(aes(x=value, color=runif_samplesize)) +
  geom_density() +
  theme_bw()

ecdf.functions <- apply(NullSimulatedTestStats, 1, ecdf)
ecdf.functions[[1]](1)

[1] 0.9036

Now calculate pi1

dat.split <- dat %>%
  filter(PC1 %in% RecodeIncludePCs) %>%
  filter(PC2 %in% RecodeIncludePCs) %>%
  mutate(PC1 = recode(PC1, !!!RecodeVec)) %>%
  mutate(PC2 = recode(PC2, !!!RecodeVec)) %>%
  group_by(PC1, P1, PC2) %>%
  mutate(test.stat.obs = -log10(min(trait.x.p.in.y))) %>%
  ungroup() %>%
  add_count(PC1, P1, PC2) %>%
  filter(n<=100) %>%
  group_by(PC1, PC2) %>%
  rowwise() %>%
  mutate(Pvals.For.Pi1 = 1-ecdf.functions[[n]](test.stat.obs)) %>%
  ungroup() %>%
  select(PC1, PC2, Pvals.For.Pi1) %>%
  filter(!PC1==PC2) %>%
  split(paste(.$PC1, .$PC2, sep = ";"))

dat.pi1 <- lapply(dat.split, CalculatePi1) %>%
  unlist() %>%
  data.frame(pi1=.) %>%
  rownames_to_column("PC1_PC2") %>%
  separate(PC1_PC2, into=c("PC1", "PC2"), sep=';')

pi.heatmap <- ggplot(dat.pi1, aes(x=PC1, y=PC2, fill=pi1)) +
  geom_raster() +
  geom_text(aes(label=signif(pi1*100, 2)), color="blue") +
  scale_fill_viridis(option="A", direction = -1, limits=c(0,1)) +
  coord_flip() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="Discovery QTL phenotype", y="Phenotype assessed for overlap")
pi.heatmap  

Now make the same heatmap but compare discovery FDR of 0.1 (like above), 0.05, 0.01

# Make a function to re-use

CalculatePi1Matrix <- function(dat, DiscoveryFDR=0.1){
  dat.split <- dat %>%
  filter(PC1 %in% RecodeIncludePCs) %>%
  filter(PC2 %in% RecodeIncludePCs) %>%
  filter(FDR.x <= DiscoveryFDR) %>%
  mutate(PC1 = recode(PC1, !!!RecodeVec)) %>%
  mutate(PC2 = recode(PC2, !!!RecodeVec)) %>%
  group_by(PC1, P1, PC2) %>%
  mutate(test.stat.obs = -log10(min(trait.x.p.in.y))) %>%
  ungroup() %>%
  add_count(PC1, P1, PC2) %>%
  filter(n<=100) %>%
  group_by(PC1, PC2) %>%
  rowwise() %>%
  mutate(Pvals.For.Pi1 = 1-ecdf.functions[[n]](test.stat.obs)) %>%
  ungroup() %>%
  select(PC1, PC2, Pvals.For.Pi1) %>%
  filter(!PC1==PC2) %>%
  split(paste(.$PC1, .$PC2, sep = ";"))

dat.pi1 <- lapply(dat.split, CalculatePi1) %>%
  unlist() %>%
  data.frame(pi1=.) %>%
  rownames_to_column("PC1_PC2") %>%
  separate(PC1_PC2, into=c("PC1", "PC2"), sep=';') %>%
  mutate(DiscoveryFDR = DiscoveryFDR)
return(dat.pi1)
}

dat.pi.AtThresholds <- bind_rows(
  CalculatePi1Matrix(dat, DiscoveryFDR=0.1),
  CalculatePi1Matrix(dat, DiscoveryFDR=0.05),
  CalculatePi1Matrix(dat, DiscoveryFDR=0.01),
)

pi.heatmap <- ggplot(dat.pi.AtThresholds, aes(x=PC1, y=PC2, fill=pi1)) +
  geom_raster() +
  geom_text(aes(label=signif(pi1*100, 2)), color="blue", size=2) +
  scale_fill_viridis(option="A", direction = -1, limits=c(0,1)) +
  facet_wrap(~DiscoveryFDR) +
  coord_flip() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x="Discovery QTL phenotype", y="Phenotype assessed for overlap", "Pi1 of xQTL among Discovery QTL classes")
pi.heatmap

Colocalization Rate

dat.coloc <- Sys.glob("../code/hyprcoloc/Results/ForColoc/MolColocTest*_*/results.txt.gz") %>%
  setNames(str_replace(., "../code/hyprcoloc/Results/ForColoc/MolColocTest(.*?)_(.+?)/results.txt.gz", "\\1_0.\\2")) %>%
  lapply(read_tsv) %>%
  bind_rows(.id="Threshold")

dat.coloc.tidy <- dat.coloc %>%
  unite(Locus, GeneLocus, Threshold, sep = ":") %>%
  left_join(., ., by = "Locus") %>% 
    filter(Trait.x != Trait.y) %>% 
    separate(Trait.x, into = c("PC1","P1"), sep = "[, ;]", remove=F) %>% 
    separate(Trait.y, into = c("PC2","P2"), sep = "[, ;]", remove=F) %>% 
    rowwise() %>%
    mutate(PC_ClassPair = paste(PC1, PC2)) %>%
    ungroup() %>%
    # pull(PC_ClassPair) %>% unique()
    mutate(IsColocalizedPair = HyprcolocIteration.x == HyprcolocIteration.y) %>%
    replace_na(list(IsColocalizedPair = FALSE))

NumQTLs <- dat.coloc.tidy %>%
  filter(PC1 %in% RecodeIncludePCs) %>%
  mutate(PC1 = recode(PC1, !!!RecodeVec)) %>%
  mutate(PC2 = recode(PC2, !!!RecodeVec)) %>%  separate(Locus, into=c("Locus", "Threshold"), sep = ":") %>%
  distinct(Locus, Threshold, PC1, P1) %>%
  count(Threshold, PC1)

coloc.rate.dat <- dat.coloc.tidy %>%
  filter(PC1 %in% RecodeIncludePCs) %>%
  filter(PC2 %in% RecodeIncludePCs) %>%
  mutate(PC1 = recode(PC1, !!!RecodeVec)) %>%
  mutate(PC2 = recode(PC2, !!!RecodeVec)) %>%
  filter(IsColocalizedPair) %>%
  separate(Locus, into=c("Locus", "Threshold"), sep = ":") %>%
  mutate(PC2 = as.factor(PC2)) %>%
  count(Locus, Threshold, PC1, P1, PC2, .drop=F) %>%
  mutate(ContainsAtleastOneColoc = n>0) %>%
  group_by(PC1, PC2, Threshold) %>%
  summarise(SumPC1sWithAtLeast1PC2Coloc = sum(ContainsAtleastOneColoc)) %>%
  left_join(NumQTLs, by=c("Threshold", "PC1")) %>%
  mutate(FractionColocs = SumPC1sWithAtLeast1PC2Coloc/n)

coloc.rate.dat %>%
  mutate(Threshold = str_remove(Threshold, "^_")) %>%
  filter(Threshold %in% c("0.001", "0.005", "0.01")) %>%
ggplot(aes(x=PC1, y=PC2, fill=FractionColocs)) +
  geom_raster() +
  geom_text(aes(label=signif(FractionColocs*100, 2)), color="blue", size=2) +
  scale_fill_viridis(option="A", direction = -1, limits=c(0,1)) +
  facet_wrap(~Threshold) +
  coord_flip() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  facet_wrap(~Threshold) +
  labs(x="Discovery QTL phenotype", y="Phenotype assessed for colocalization", "Percent Discovery QTLs with at least one colocalizing xQTL")

### ColocalizationRateAmongAttempted

coloc.rate.AmongAttempted.dat <- dat.coloc.tidy %>%
  filter(PC1 %in% RecodeIncludePCs) %>%
  filter(PC2 %in% RecodeIncludePCs) %>%
  mutate(PC1 = recode(PC1, !!!RecodeVec)) %>%
  mutate(PC2 = recode(PC2, !!!RecodeVec)) %>%
  separate(Locus, into=c("Locus", "Threshold"), sep = ":") %>%
  group_by(PC1, PC2, Threshold) %>%
  summarise(FractionColocs = sum(IsColocalizedPair)/n())

coloc.rate.AmongAttempted.dat %>%
  mutate(Threshold = str_remove(Threshold, "^_")) %>%
  filter(Threshold %in% c("0.001", "0.005", "0.01")) %>%
ggplot(aes(x=PC1, y=PC2, fill=FractionColocs)) +
  geom_raster() +
  geom_text(aes(label=signif(FractionColocs*100, 2)), color="blue", size=2) +
  scale_fill_viridis(option="A", direction = -1, limits=c(0,1)) +
  facet_wrap(~Threshold) +
  coord_flip() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  facet_wrap(~Threshold) +
  labs(x="Phenotype_x", y="Phenotype_y", title="Percent succesful colocalizations among attempted")