Thesis MGRB Data Comparison Script

ALL IMPORTED FILES MUST BE IN THE SAME DIRECTORY (‘MGRB’) AS THIS SCRIPT

Load required packages

library("readbulk")
library("tidyverse")

Import files with sample cohort and relevant MGRB data for individual candidate genes (located in sub-directory ‘Genes’, downloaded from https://sgc.garvan.org.au/), removing other genes

ViP_list <- read.delim("ViP Samples & Candidate Genes.txt", header=TRUE, stringsAsFactors=FALSE, sep="\t", comment.char="") %>% 
  arrange(Sample)
genes <- ViP_list$SYMBOL

MGRB <- read_bulk(directory="Genes")
MGRB_genes <- filter(MGRB, gene%in%genes) %>% 
  filter(!gene%in%"ENTPD4") %>% 
  mutate_at(vars(gnomadAF),funs(replace(., is.na(.), 0)))

Remove non-LOF and common (GnomAD AF > 0.005) variants

MGRB_genes_LOF <- filter(MGRB_genes,consequences%in%c(
  "frameshift_variant",
  "stop_gained",
  "splice_donor_variant",
  "splice_acceptor_variant",
  "start_lost",
  "frameshift_variant,splice_region_variant"))
MGRB_genes_LOF_rare <- filter(MGRB_genes_LOF, gnomadAF<=0.005)

Calculate summary stats for MGRB LoF candidate gene variants, and output results

MGRB_LOF_stats <-group_by(MGRB_genes_LOF,gene) %>%
  select(Allele.Count,Allele.Frequency) %>% 
  mutate(Total_MGRB_Allele_Count=sum(Allele.Count)) %>% 
  mutate(Total_MGRB_Gene_Freq=sum(Allele.Frequency)) %>% 
  distinct(gene,Total_MGRB_Allele_Count,Total_MGRB_Gene_Freq) %>% 
  rename(SYMBOL=gene)
MGRB_LOF_rare_stats <-group_by(MGRB_genes_LOF_rare,gene) %>%
  select(Allele.Count,Allele.Frequency) %>% 
  mutate(Total_MGRB_Allele_Count_0.005=sum(Allele.Count)) %>% 
  mutate(Total_MGRB_Gene_Freq_0.005=sum(Allele.Frequency)) %>% 
  distinct(gene,Total_MGRB_Allele_Count_0.005,Total_MGRB_Gene_Freq_0.005) %>% 
  rename(SYMBOL=gene)

write_excel_csv(MGRB_LOF_stats,path="MGRB_LOF_stats.csv", na=".",append=FALSE,col_names=TRUE)
write_excel_csv(MGRB_LOF_rare_stats,path="MGRB_LOF_rare_stats.csv", na=".",append=FALSE,col_names=TRUE)

Import candidate gene summary stats for case cohort (see Nature Comms paper- https://www.nature.com/articles/s41467-020-15461-z, and Supp Data 2), and merge with MGRB summary data

ViP_Candidate_Genes_Stats <- read.delim("ViP Candidate Genes Stats.txt", header=TRUE, stringsAsFactors=FALSE, sep="\t", comment.char="")
ViP_MGRB_Stats_Combined <-  left_join(ViP_Candidate_Genes_Stats, MGRB_LOF_stats,by="SYMBOL",copy=FALSE) %>% 
  mutate_at(vars(Total_MGRB_Allele_Count),funs(replace(., is.na(.), 0))) %>% 
  mutate_at(vars(Total_MGRB_Gene_Freq),funs(replace(., is.na(.), 0)))
ViP_MGRB_Stats_Combined <- left_join(ViP_MGRB_Stats_Combined, MGRB_LOF_rare_stats,by="SYMBOL",copy=FALSE) %>% 
  mutate_at(vars(Total_MGRB_Allele_Count_0.005),funs(replace(., is.na(.), 0))) %>% 
  mutate_at(vars(Total_MGRB_Gene_Freq_0.005),funs(replace(., is.na(.), 0))) %>% 
  mutate(Total_MGRB_Alleles=5690)

Calculate Fisher’s test values, and output results

ViP_MGRB_Stats_Combined_ratios_fisherresults <- 
  mutate(ViP_MGRB_Stats_Combined, 
         "Sample_Gene_LOF_Freq_Ratio_MGRB"=Sample_Gene_Freq/Total_MGRB_Gene_Freq,
         "Sample_Gene_LOF_Freq_Ratio_MGRB_0.005"=Sample_Gene_Freq/Total_MGRB_Gene_Freq_0.005,
         "Ratio_Difference"=Sample_Gene_LOF_Freq_Ratio_MGRB_0.005-Sample_Gene_LOF_Freq_Ratio_MGRB)

fisherresults_MGRB_0.005 <- apply(ViP_MGRB_Stats_Combined_ratios_fisherresults[,c(16,29,18,31)],1,function(x) 
  {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})
fisherresults_MGRBvsGnomADNFE_0.005 <- apply(ViP_MGRB_Stats_Combined_ratios_fisherresults[,c(29,17,31,19)],1,function(x) 
{fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

ViP_MGRB_Stats_Combined_ratios_fisherresults$P_value_Fisher_MGRB_0.005 = sapply(fisherresults_MGRB_0.005,function(x) round(x$p.value,12))
ViP_MGRB_Stats_Combined_ratios_fisherresults$OR_MGRB_0.005 = sapply(fisherresults_MGRB_0.005,function(x) round(x$estimate,4))
ViP_MGRB_Stats_Combined_ratios_fisherresults$`95%CI_MGRB_0.005` = sapply(fisherresults_MGRB_0.005,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
ViP_MGRB_Stats_Combined_ratios_fisherresults$P_value_Fisher_MGRBvsGnomADNFE_0.005 = sapply(fisherresults_MGRBvsGnomADNFE_0.005,function(x) round(x$p.value,12))
ViP_MGRB_Stats_Combined_ratios_fisherresults$OR_MGRBvsGnomADNFE_0.005 = sapply(fisherresults_MGRBvsGnomADNFE_0.005,function(x) round(x$estimate,4))
ViP_MGRB_Stats_Combined_ratios_fisherresults$`95%CI_MGRBvsGnomADNFE_0.005` = sapply(fisherresults_MGRBvsGnomADNFE_0.005,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

write_excel_csv(ViP_MGRB_Stats_Combined_ratios_fisherresults,path="ViP_MGRB_Stats_Combined_ratios_fisherresults.csv", na=".",append=FALSE,col_names=TRUE)