Filtering and ranking R script for ‘Exome sequencing of familial high-grade serous ovarian carcinoma reveals heterogeneity for rare candidate susceptibility genes’ (Subramanian et al., Nature Communications, 2020)

ALL IMPORTED FILES MUST BE IN THE SAME DIRECTORY AS THIS SCRIPT

Load tidyverse package

library("tidyverse")

Import raw vcf data from “masterfile” (vcf file output from alignment, variant caller and annotation pipeline- available from authors on reasonable request) and exclude variants with GnomAD non-cancer AF > 0.005

masterfile_eoc <- read.delim("masterfile.tsv", header=TRUE, row.names=NULL, na.strings = ".", stringsAsFactors=FALSE) %>% 
  dplyr::rename("Sample"="X.Sample") %>% 
  filter(GnomAD_v2.1_non_cancer_AF<=0.005)

Append patient pathology data (available from authors on reasonable request)

OvCa_Path_Data <- read.delim("OvCa_Path_Data.txt", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>%
  dplyr::rename("Sample"=Exome.ID)
masterfile_eoc2 <- left_join(masterfile_eoc,OvCa_Path_Data,by="Sample",copy=FALSE)
rm(OvCa_Path_Data)

Exclude low-quality variants and GnomAD RF/InbreedingCoeff-flagged variants

masterfile_eoc_goodQ <- filter(masterfile_eoc2, (QUAL>=30)&
                                 (Identified!="FilteredInAll")&
                                 (Sample.PMCDP>=10)&
                                 (Sample.PMCFREQ>=0.25)) %>% 
                                 filter((!str_detect(GnomAD_v2.1_FILTER_exome,"InbreedingCoeff")&
                                           !str_detect(GnomAD_v2.1_FILTER_exome,"RF"))|
                                       is.na(GnomAD_v2.1_FILTER_exome)) %>% 
  filter((!str_detect(GnomAD_v2.1_FILTER_genome,"InbreedingCoeff")&
            !str_detect(GnomAD_v2.1_FILTER_genome,"RF"))|
           is.na(GnomAD_v2.1_FILTER_genome))

Exclude other mutation +ve samples (MLH1/MSH2/MSH6/PMS2, TP53, RAD51C/D, BRIP1)

masterfile_eoc_goodQ2 <- filter(masterfile_eoc_goodQ,
                                (Sample!="PUB-P1RUC")&
                                  (Sample!="PUB-06NJ4")&
                                  (Sample!="PUB-4WIFF")&
                                  (Sample!="PUB-JYFJG")&
                                  (Sample!="PUB-KJM27")&
                                  (Sample!="PUB-RT7RU"))

Exclude VEP MODERATE impact variants

masterfile_eoc_goodQ2_HIGH <- filter(masterfile_eoc_goodQ2, IMPACT=="HIGH")

Separate Ensembl canonical variants

masterfile_eoc_goodQ2_HIGH_ENSTcanonical <- filter(masterfile_eoc_goodQ2_HIGH,CANONICAL=="YES") 

Sample variant counts and AFs

variants_heterozygous <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical %>% 
  select(HGVSc,Sample.GT) %>% 
  group_by(HGVSc,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT%in%c("'0/1","'1/0"))) %>% 
  rename(alleles = "n()") %>% 
  group_by(HGVSc) %>% 
  summarise(sum(alleles))
variants_homozygous <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical %>% 
  select(HGVSc,Sample.GT) %>% 
  group_by(HGVSc,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT!="'0/1")&(Sample.GT!="'1/0")) %>% 
  mutate(alleles = `n()`*2) %>% 
  select(-`n()`) %>% 
  group_by(HGVSc) %>% 
  summarise(sum(alleles))
variants <- full_join(variants_heterozygous,variants_homozygous,by="HGVSc",copy=FALSE,suffix=c(".x",".y")) %>%
  mutate_all(funs(replace(., is.na(.), 0))) %>% 
  mutate(Total_Allele_Count=`sum(alleles).x`+`sum(alleles).y`) %>% 
  mutate(Sample_AF=Total_Allele_Count/(n_distinct(masterfile_eoc_goodQ2$Sample)*2))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_2 <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical,variants[,c(1,4:5)],by="HGVSc",copy=FALSE)

Exclude variants with sample AF > 0.01̛

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_2,Sample_AF < 0.01)

Sample gene counts and frequencies

genes_oneVarAllele <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 %>%
  select(SYMBOL,Gene,Sample.GT) %>% 
  group_by(Gene,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT%in%c("'0/1","'1/0"))) %>% 
  rename(gene_Var = "n()") %>% 
  group_by(Gene) %>% 
  summarise(sum(gene_Var))
genes_twoVarAllele <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 %>%
  select(SYMBOL,Gene,Sample.GT) %>% 
  group_by(Gene,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT!="'0/1")&(Sample.GT!="'1/0")) %>% 
  mutate(gene_Var = `n()`*2) %>% 
  select(-`n()`) %>% 
  group_by(Gene) %>% 
  summarise(sum(gene_Var))
genes <- full_join(genes_oneVarAllele,genes_twoVarAllele,by="Gene",copy=FALSE,suffix=c(".x",".y")) %>% 
  mutate_all(funs(replace(., is.na(.), 0))) %>% 
  mutate(Total_Gene_Count=`sum(gene_Var).x`+`sum(gene_Var).y`) %>% 
  mutate(Sample_Gene_Freq=Total_Gene_Count/(n_distinct(masterfile_eoc_goodQ2$Sample)*2))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2 <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01, genes[,c(1,4:5)],by="Gene",copy=FALSE)

Add GnomAD aggregated gene-level variant frequency data (available from authors on reasonable request)

GnomAD_stats_LOF_VEP <- read.delim("GnomAD_stats_LOF_VEP.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE)

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2, GnomAD_stats_LOF_VEP[,c(2,30:41)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats2 <- mutate_at(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats,vars(starts_with("FILTER_")),funs(replace(., is.na(.), 0))) %>% 
  mutate_if(grepl("popmax$", names(.)),funs(ifelse(. == "NA", 0, as.numeric(.)))) %>% 
  mutate_at(vars(ends_with("popmax")),funs(replace(., is.na(.), 0)))

Calculate Sample/GnomAD gene frequency ratios for total and NFE GnomAD figures, including and excluding variants with AF > 0.005

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats2, 
         "Sample_Gene_LOF_Freq_Ratio_GnomAD"=Sample_Gene_Freq/FILTER_RF_LOF_HIGHIMPACT_AF_1.0_ADJ,
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_0.005"=Sample_Gene_Freq/FILTER_RF_LOF_HIGHIMPACT_AF_0.005_ADJ,
         "Ratio_Difference"=Sample_Gene_LOF_Freq_Ratio_GnomAD_0.005-Sample_Gene_LOF_Freq_Ratio_GnomAD, 
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE"=Sample_Gene_Freq/FILTER_RF_LOF_HIGHIMPACT_AF_1.0_NFE_ADJ,
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.005"=Sample_Gene_Freq/FILTER_RF_LOF_HIGHIMPACT_AF_0.005_NFE_ADJ,
         "Ratio_Difference_NFE"=Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.005-Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE)

Calculate Sample/GnomAD allele frequency ratios

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2 <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios, "Sample_Gene_LOF_AF_Ratio"=Sample_AF/GnomAD_v2.1_non_cancer_AF, 
"Sample_Gene_LOF_AF_Ratio_NFE"=Sample_AF/GnomAD_v2.1_non_cancer_AF_nfe)

Create data frames with Agilent SureSelect whole exome genes (all ENST, protein-coding only and non-protein-coding- files available from authors on reasonable request) for waterfall plot

ensembl_biotypes <- read.delim("ensembl_biotypes.tsv", stringsAsFactors=FALSE) %>% 
  select(ensembl_transcript_id,transcript_biotype) %>% 
  dplyr::rename("Feature"="ensembl_transcript_id","BIOTYPE"="transcript_biotype")
GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all <- read.delim("GnomAD.gene.stats.v2.RF.only.agilent.v6.only.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  filter(CANONICAL=="YES") %>% 
  dplyr::rename("SYMBOL"="X.Symbol") %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE)
GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding <- read.delim("GnomAD.gene.stats.v2.RF.only.agilent.v6.only.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  filter(CANONICAL=="YES") %>% 
  dplyr::rename("SYMBOL"="X.Symbol") %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE) %>% 
  filter(BIOTYPE%in%c("protein_coding"))
GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding <- read.delim("GnomAD.gene.stats.v2.RF.only.agilent.v6.only.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  filter(CANONICAL=="YES") %>% 
  dplyr::rename("SYMBOL"="X.Symbol") %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE) %>% 
  filter(!BIOTYPE%in%c("protein_coding"))

allGenes_AgilentSSv6_list <- tibble("SYMBOL"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all$SYMBOL,"Gene"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all$Gene,"Feature"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all$Feature) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)
allGenes_AgilentSSv6_list_protein_coding_only <- tibble("SYMBOL"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding$SYMBOL,"Gene"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding$Gene,"Feature"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding$Feature) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)
allGenes_AgilentSSv6_list_non_protein_coding_only <- tibble("SYMBOL"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding$SYMBOL,"Gene"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding$Gene,"Feature"=GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding$Feature) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)

allGenes_masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_list <- tibble("SYMBOL"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$SYMBOL,"Gene"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$Gene,"BIOTYPE"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)
allGenes_masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_list <- tibble("SYMBOL"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$SYMBOL,"Gene"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$Gene,"BIOTYPE"=masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  arrange(SYMBOL)

waterfall_sampleAF0.01_allGenes <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2 %>% 
  mutate(Total_Case_Alleles=(n_distinct(masterfile_eoc_goodQ2$Sample)*2)) %>% 
  select(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) %>% 
  distinct(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) 
waterfall_sampleAF0.01_allGenes2 <- bind_rows(waterfall_sampleAF0.01_allGenes, anti_join(allGenes_AgilentSSv6_list,waterfall_sampleAF0.01_allGenes,by="Gene")) %>% 
  mutate_at(vars("Total_Gene_Count"),funs(replace(., is.na(.), 0))) %>% 
  mutate_at(vars("Total_Case_Alleles"),funs(replace(., is.na(.), n_distinct(masterfile_eoc_goodQ2$Sample)*2)))
waterfall_sampleAF0.01_allGenes3 <- right_join(
  waterfall_sampleAF0.01_allGenes2,
  GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all,
  by="Gene",copy=FALSE) %>% 
  select("SYMBOL.x","Gene","Total_Gene_Count","Total_Case_Alleles","FILTER_RF_LOF_HIGHIMPACT_AC_0.005","FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE","MAX_AN","MAX_AN_NFE","FILTER_RF_LOF_HIGHIMPACT_AF_0.005_ADJ","FILTER_RF_LOF_HIGHIMPACT_AF_0.005_NFE_ADJ") %>% 
  dplyr::rename("SYMBOL"="SYMBOL.x") %>%
  mutate(FILTER_RF_LOF_HIGHIMPACT_AC_0.005_ADJ=round((MAX_AN*FILTER_RF_LOF_HIGHIMPACT_AF_0.005_ADJ))) %>% 
  mutate(FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE_ADJ=round((MAX_AN_NFE*FILTER_RF_LOF_HIGHIMPACT_AF_0.005_NFE_ADJ))) %>% 
  distinct(.keep_all=TRUE)

waterfall_sampleAF0.01_biotype_allGenes <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2 %>% 
  mutate(Total_Case_Alleles=(n_distinct(masterfile_eoc_goodQ2$Sample)*2)) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  select(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) %>% 
  distinct(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) 
waterfall_sampleAF0.01_biotype_allGenes2 <- bind_rows(waterfall_sampleAF0.01_biotype_allGenes, anti_join(allGenes_AgilentSSv6_list_protein_coding_only,waterfall_sampleAF0.01_biotype_allGenes,by="Gene")) %>% 
  mutate_at(vars("Total_Gene_Count"),funs(replace(., is.na(.), 0))) %>% 
  mutate_at(vars("Total_Case_Alleles"),funs(replace(., is.na(.), n_distinct(masterfile_eoc_goodQ2$Sample)*2)))
waterfall_sampleAF0.01_biotype_allGenes3 <- right_join(
  waterfall_sampleAF0.01_biotype_allGenes2,
  GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding,
  by="Gene",copy=FALSE) %>% 
  select("SYMBOL.x","Gene","Total_Gene_Count","Total_Case_Alleles","FILTER_RF_LOF_HIGHIMPACT_AC_0.005","FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE","MAX_AN","MAX_AN_NFE","FILTER_RF_LOF_HIGHIMPACT_AF_0.005_ADJ","FILTER_RF_LOF_HIGHIMPACT_AF_0.005_NFE_ADJ") %>% 
  dplyr::rename("SYMBOL"="SYMBOL.x") %>% 
  mutate(FILTER_RF_LOF_HIGHIMPACT_AC_0.005_ADJ=round((MAX_AN*FILTER_RF_LOF_HIGHIMPACT_AF_0.005_ADJ))) %>% 
  mutate(FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE_ADJ=round((MAX_AN_NFE*FILTER_RF_LOF_HIGHIMPACT_AF_0.005_NFE_ADJ))) %>% 
  distinct(.keep_all=TRUE)

Calculate Fisher’s test values

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios2,Total_Case_Alleles=(n_distinct(masterfile_eoc_goodQ2$Sample)*2)) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, na.rm=TRUE))) %>% 
  mutate(MAX_AN_NFE=replace(.$MAX_AN_NFE, is.na(.$MAX_AN_NFE), max(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_ENST <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(221,229,243,244)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005 = sapply(fisherresults_ENST,function(x) round(x$p.value,12))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005 = sapply(fisherresults_ENST,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$`95%CI_0.005` = sapply(fisherresults_ENST,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults2_ENST <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(221,232,243,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005_NFE = sapply(fisherresults2_ENST,function(x) round(x$p.value,12))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005_NFE = sapply(fisherresults2_ENST,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$`95%CI_0.005_NFE` = sapply(fisherresults2_ENST,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

Calculate Fisher’s test values (waterfall plot)

fisherresults_allGenes_waterfall <- apply(waterfall_sampleAF0.01_allGenes3[,c(3,5,4,7)],1,function(x) {fisher.test(rbind(x[1:2],c(abs(x[3]-x[1]),abs(x[4]-x[2]))))})

waterfall_sampleAF0.01_allGenes3$P_value_Fisher_0.005 = sapply(fisherresults_allGenes_waterfall,function(x) round(x$p.value,15))
waterfall_sampleAF0.01_allGenes3$OR_0.005 = sapply(fisherresults_allGenes_waterfall,function(x) round(x$estimate,15))
waterfall_sampleAF0.01_allGenes3$`95%CI_0.005` = sapply(fisherresults_allGenes_waterfall,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_NFE_waterfall <- apply(waterfall_sampleAF0.01_allGenes3[,c(3,6,4,8)],1,function(x) {fisher.test(rbind(x[1:2],c(abs(x[3]-x[1]),abs(x[4]-x[2]))))})

waterfall_sampleAF0.01_allGenes3$P_value_Fisher_0.005_NFE = sapply(fisherresults_NFE_waterfall,function(x) round(x$p.value,15))
waterfall_sampleAF0.01_allGenes3$OR_0.005_NFE = sapply(fisherresults_NFE_waterfall,function(x) round(x$estimate,15))
waterfall_sampleAF0.01_allGenes3$`95%CI_0.005_NFE` = sapply(fisherresults_NFE_waterfall,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_biotype_allGenes_waterfall <- apply(waterfall_sampleAF0.01_biotype_allGenes3[,c(3,5,4,7)],1,function(x) {fisher.test(rbind(x[1:2],c(abs(x[3]-x[1]),abs(x[4]-x[2]))))})

waterfall_sampleAF0.01_biotype_allGenes3$P_value_Fisher_0.005 = sapply(fisherresults_biotype_allGenes_waterfall,function(x) round(x$p.value,15))
waterfall_sampleAF0.01_biotype_allGenes3$OR_0.005 = sapply(fisherresults_biotype_allGenes_waterfall,function(x) round(x$estimate,15))
waterfall_sampleAF0.01_biotype_allGenes3$`95%CI_0.005` = sapply(fisherresults_biotype_allGenes_waterfall,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_biotype_NFE_waterfall <- apply(waterfall_sampleAF0.01_biotype_allGenes3[,c(3,6,4,8)],1,function(x) {fisher.test(rbind(x[1:2],c(abs(x[3]-x[1]),abs(x[4]-x[2]))))})

waterfall_sampleAF0.01_biotype_allGenes3$P_value_Fisher_0.005_NFE = sapply(fisherresults_biotype_NFE_waterfall,function(x) round(x$p.value,15))
waterfall_sampleAF0.01_biotype_allGenes3$OR_0.005_NFE = sapply(fisherresults_biotype_NFE_waterfall,function(x) round(x$estimate,15))
waterfall_sampleAF0.01_biotype_allGenes3$`95%CI_0.005_NFE` = sapply(fisherresults_biotype_NFE_waterfall,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

Annotate variants with filtering AF > maximum credible population AF in GnomAD for hereditary OvCa (0.000181), calculated using alleleFrequencyApp (http://cardiodb.org/allelefrequencyapp/)

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2, FilteringAF_95_Exceeds_MaxCredPopAF=((GnomAD_v2.1_non_cancer_Filtering_AF_95)>0.000181)) %>% 
  mutate(FilteringAF_95_NFE_Exceeds_MaxCredPopAF=((GnomAD_v2.1_non_cancer_Filtering_AF_95_nfe)>0.000181))

Annotate variants that PASS/FAIL the LOFTEE 50bp rule

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4 <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3,LoF_50_BP_RULE = ifelse(grepl("50_BP_RULE:PASS", LoF_info), TRUE, ifelse(grepl("50_BP_RULE:FAIL", LoF_info), FALSE, NA)))

Import NMD depleted gene list (available from authors on reasonable request) and annotate variants accordingly

NMD_depleted_gene_list <- read.csv("NMD_depleted_gene_list.csv", stringsAsFactors=FALSE) %>% 
  rename("SYMBOL" = gene) %>% 
  rename("Feature" = txnames) %>% 
  rename("NMD_predictor_rank" = min.rank)

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults5 <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4, NMD_depleted_gene_list[,c(2:3)],by="Feature",copy=FALSE) 

Output data (pre-prioritisation) and ranked p-value list used for Benjamini-Hochberg adjustment

write_excel_csv(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults5,path="masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01.csv",na=".",append=FALSE,col_names=TRUE)

select(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults5,SYMBOL,Gene,P_value_Fisher_0.005_NFE) %>% 
  arrange(P_value_Fisher_0.005_NFE) %>% 
  unique() %>% 
  write_tsv(path="HIGH_ENSTcanonical_BH_p_values.tsv")

Calculate total LoF variants in sample vs total LoF variants in GnomAD non-cancer (all/NFE-only) for all genes and protein-coding-only genes, and use for chi-squared test

a <- sum(waterfall_sampleAF0.01_allGenes3$Total_Gene_Count) %>% as.numeric()
b <- n_distinct(masterfile_eoc_goodQ2$Sample)*2*(nrow(allGenes_AgilentSSv6_list)) %>% as.numeric()
b1 <- b-a %>% as.numeric()
c <- sum(waterfall_sampleAF0.01_allGenes3$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
d <- (max(waterfall_sampleAF0.01_allGenes3$MAX_AN))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
d1 <- d-c %>% as.numeric()
e <- sum(waterfall_sampleAF0.01_allGenes3$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
f <- (max(waterfall_sampleAF0.01_allGenes3$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
f1 <- f-e %>% as.numeric()

g <- sum(waterfall_sampleAF0.01_biotype_allGenes3$Total_Gene_Count) %>% as.numeric()
h <- n_distinct(masterfile_eoc_goodQ2$Sample)*2*(nrow(allGenes_AgilentSSv6_list_protein_coding_only)) %>% as.numeric()
h1 <- h-g %>% as.numeric()
i <- sum(waterfall_sampleAF0.01_biotype_allGenes3$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
j <- (max(waterfall_sampleAF0.01_biotype_allGenes3$MAX_AN))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
j1 <- j-i %>% as.numeric()
k <- sum(waterfall_sampleAF0.01_biotype_allGenes3$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
l <- (max(waterfall_sampleAF0.01_biotype_allGenes3$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
l1 <- l-k %>% as.numeric()

chiX2_0.005 <- matrix(c(a,b1,c,d1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE <- matrix(c(a,b1,e,f1),nrow=2,byrow=TRUE)
chisq.test(chiX2_0.005,correct=FALSE)
chisq.test(chiX2_0.005_NFE,correct=FALSE)

chiX2_0.005_biotype <- matrix(c(g,h1,i,j1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_biotype <- matrix(c(g,h1,k,l1),nrow=2,byrow=TRUE)
chisq.test(chiX2_0.005_biotype,correct=FALSE)
chisq.test(chiX2_0.005_NFE_biotype,correct=FALSE)

Remove variants with non-protein-coding biotypes

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults5,BIOTYPE%in%c("protein_coding"))

Remove variants with GnomAD PopMax AF > 0.005

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005 <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype,
GnomAD_v2.1_non_cancer_AF_popmax<=0.005)

Remove variants with ratio < 3

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3 <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005,Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.005>3) 

Remove variants with gene counts < 3

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_geneCount3 <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3,
Total_Gene_Count>2)

Output data (post-prioritisation)

write_excel_csv(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_geneCount3,path="masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_geneCount3.csv",na=".",append=FALSE,col_names=TRUE)

Output data for log p-value waterfall plots (all genes and protein-coding-only genes)

waterfall_sampleAF0.01_allGenes4 <- mutate(waterfall_sampleAF0.01_allGenes3,"Log10_P-value_Fisher_0.005_NFE" = log10(P_value_Fisher_0.005_NFE)) %>% 
  arrange(desc(OR_0.005_NFE)) %>% 
  mutate("Plotted_Log10_P-value_Fisher_0.005_NFE"= 
  if_else (OR_0.005_NFE>=1, `Log10_P-value_Fisher_0.005_NFE`*-1, `Log10_P-value_Fisher_0.005_NFE`, missing=NULL)) %>% 
  arrange(desc(`Plotted_Log10_P-value_Fisher_0.005_NFE`),desc(OR_0.005_NFE)) %>% 
  write_excel_csv(path="waterfall_sampleAF0.01_allGenes_v6.csv",na=".",append=FALSE,col_names=TRUE)

waterfall_sampleAF0.01_biotype_allGenes4 <- mutate(waterfall_sampleAF0.01_biotype_allGenes3,"Log10_P-value_Fisher_0.005_NFE" = log10(P_value_Fisher_0.005_NFE)) %>% 
  arrange(desc(OR_0.005_NFE)) %>% 
  mutate("Plotted_Log10_P-value_Fisher_0.005_NFE"= 
  if_else (OR_0.005_NFE>=1, `Log10_P-value_Fisher_0.005_NFE`*-1, `Log10_P-value_Fisher_0.005_NFE`, missing=NULL)) %>% 
  arrange(desc(`Plotted_Log10_P-value_Fisher_0.005_NFE`),desc(OR_0.005_NFE)) %>% 
  write_excel_csv(path="waterfall_sampleAF0.01_biotype_allGenes_v6.csv",na=".",append=FALSE,col_names=TRUE)

Remove redundant c.2208T>G stopgain BLM variant from PUB-1G5DH (in cis with c.2206dupT variant in this individual) and recalculate gene counts/frequencies and ratios for DNA repair gene analysis

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_a <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01, HGVSc!="ENST00000355112.3:c.2208T>G")

genes_oneVarAllele <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_a %>% 
  select(SYMBOL,Sample.GT) %>% 
  group_by(SYMBOL,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT%in%c("'0/1","'1/0"))) %>% 
  rename(gene_Var = "n()") %>% 
  group_by(SYMBOL) %>% 
  summarise(sum(gene_Var))
genes_twoVarAllele <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_a %>% 
  select(SYMBOL,Sample.GT) %>% 
  group_by(SYMBOL,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT!="'0/1")&(Sample.GT!="'1/0")) %>% 
  mutate(gene_Var = `n()`*2) %>% 
  select(-`n()`) %>% 
  group_by(SYMBOL) %>% 
  summarise(sum(gene_Var))
genes <- full_join(genes_oneVarAllele,genes_twoVarAllele,by="SYMBOL",copy=FALSE,suffix=c(".x",".y")) %>% 
  mutate_all(funs(replace(., is.na(.), 0))) %>% 
  mutate(Total_Gene_Count=`sum(gene_Var).x`+`sum(gene_Var).y`) %>% 
  mutate(Sample_Gene_Freq=Total_Gene_Count/(n_distinct(masterfile_eoc_goodQ2$Sample)*2))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2a <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_a, genes[,c(1,4:5)],by="SYMBOL",copy=FALSE)

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_a <- left_join(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2a, GnomAD_stats_LOF_VEP[,c(2,30:41)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a <- mutate_at(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_a,vars(starts_with("FILTER_")),funs(replace(., is.na(.), 0))) %>% 
  mutate_if(grepl("popmax$", names(.)),funs(ifelse(. == "NA", 0, as.numeric(.)))) %>% 
  mutate_at(vars(ends_with("popmax")),funs(replace(., is.na(.), 0)))

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_a <- 
  mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a,
         "Sample_Gene_LOF_Freq_Ratio_GnomAD"=Sample_Gene_Freq/masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a[,225],
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_0.005"=Sample_Gene_Freq/masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a[,231],
         "Ratio_Difference"=Sample_Gene_LOF_Freq_Ratio_GnomAD_0.005-Sample_Gene_LOF_Freq_Ratio_GnomAD, 
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE"=Sample_Gene_Freq/masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a[,228],
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.005"=Sample_Gene_Freq/masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_2a[,234],
         "Ratio_Difference_NFE"=Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.005-Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE) %>% 
  mutate(
    "Sample_Gene_LOF_AF_Ratio"=Sample_AF/GnomAD_v2.1_non_cancer_AF, 
    "Sample_Gene_LOF_AF_Ratio_NFE"=Sample_AF/GnomAD_v2.1_non_cancer_AF_nfe)

Attach other GnomAD data for DNA repair genes analysis

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_a,Total_Case_Alleles=(n_distinct(masterfile_eoc_goodQ2$Sample)*2)) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_3a <- mutate(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a, MAX_AN=replace(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a$MAX_AN, is.na(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a$MAX_AN), max(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a$MAX_AN, na.rm=TRUE))) %>% 
  mutate(MAX_AN_NFE=replace(.$MAX_AN_NFE, is.na(.$MAX_AN_NFE), max(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_2a$MAX_AN_NFE, na.rm=TRUE)))

Extract DNA repair gene variants from sample data by repair pathway, and combine into one table

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_3a,
GnomAD_v2.1_non_cancer_AF_popmax<=0.005)
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_BER <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","XRCC1","PCNA","PARP1","PARP2","APEX1","FEN1","POLE","POLD1","LIG3","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","LIG1"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_MMR <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","POLD1"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_NER <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("RPA1","RPA2","RPA3","RPA4","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC1","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_HR <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","PALB2","MUS81","EME1","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_FA <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_NHEJ <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DR <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("MGMT","ALKBH2","ALKBH3"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_ALLrepair <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("XRCC1","PARP1","PARP2","APEX1","FEN1","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19","ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1","ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1","XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1","MGMT","ALKBH2","ALKBH3"))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_OTHER <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_a,SYMBOL%in%c("PAXIP1","BLM","KMT2C","CRIP1","CDK12","BAP1","BARD1","WRN","BUB1","CENPE","ZW10","TTK","KNTC1","AURKB","POLB","POLH","POLQ","TDP1","TDP2","NUDT1","DUT","RRM2B","POLG","REV3L","MAD2L2","REV1","POLI","POLK","POLL","POLM","POLN","TREX1","TREX2","APTX","SPO11","ENDOV","UBE2A","UBE2B","RAD18","SHPRH","HLTF","RNF168","SPRTN","RNF8","RNF4","UBE2V2","UBE2N","H2AFX","CHAF1A","SETMAR","RECQL4","MPLKIP","DCLRE1A","DCLRE1B","PRPF19","RECQL","RECQL5","HELQ","RDM1","NABP2","ATRIP","MDC1","RAD1","RAD9A","HUS1","RAD17","TP53","TP53BP1","TOPBP1","CLK2","PER1"))
DNArepair_other_list <- tibble("SYMBOL"=c("XRCC1","PARP1","PARP2","APEX1","FEN1","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19","ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1","ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86","C19orf40","MUS81","EME1","XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1","MGMT","ALKBH2","ALKBH3","PAXIP1","BLM","KMT2C","CRIP1","CDK12","BAP1","BARD1","WRN","BUB1","CENPE","ZW10","TTK","KNTC1","AURKB","POLB","POLH","POLQ","TDP1","TDP2","NUDT1","DUT","RRM2B","POLG","REV3L","MAD2L2","REV1","POLI","POLK","POLL","POLM","POLN","TREX1","TREX2","APTX","SPO11","ENDOV","UBE2A","UBE2B","RAD18","SHPRH","HLTF","RNF168","SPRTN","RNF8","RNF4","UBE2V2","UBE2N","H2AFX","CHAF1A","SETMAR","RECQL4","MPLKIP","DCLRE1A","DCLRE1B","PRPF19","RECQL","RECQL5","HELQ","RDM1","NABP2","ATRIP","MDC1","RAD1","RAD9A","HUS1","RAD17","TP53","TP53BP1","TOPBP1","CLK2","PER1"))
                                                                                                           masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other <- bind_rows(
  masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_ALLrepair,
  masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_OTHER) 

Extract DNA repair gene counts from sample data by repair pathway

genes_BER <- filter(genes,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","XRCC1","PCNA","PARP1","PARP2","APEX1","FEN1","POLE","POLD1","LIG3","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","LIG1"))
genes_MMR <- filter(genes,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","POLD1"))
genes_NER <- filter(genes,SYMBOL%in%c("RPA1","RPA2","RPA3","RPA4","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC1","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19"))
genes_HR <- filter(genes,SYMBOL%in%c("ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","PALB2","MUS81","EME1","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1"))
genes_FA <- filter(genes,SYMBOL%in%c("ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1"))
genes_NHEJ <- filter(genes,SYMBOL%in%c("XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1"))
genes_DR <- filter(genes,SYMBOL%in%c("MGMT","ALKBH2","ALKBH3"))
genes_ALLrepair <- filter(genes,SYMBOL%in%c("XRCC1","PARP1","PARP2","APEX1","FEN1","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19","ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1","ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1","XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1","MGMT","ALKBH2","ALKBH3"))
genes_OTHER <- filter(genes,SYMBOL%in%c("PAXIP1","BLM","KMT2C","CRIP1","CDK12","BAP1","BARD1","WRN","BUB1","CENPE","ZW10","TTK","KNTC1","AURKB","POLB","POLH","POLQ","TDP1","TDP2","NUDT1","DUT","RRM2B","POLG","REV3L","MAD2L2","REV1","POLI","POLK","POLL","POLM","POLN","TREX1","TREX2","APTX","SPO11","ENDOV","UBE2A","UBE2B","RAD18","SHPRH","HLTF","RNF168","SPRTN","RNF8","RNF4","UBE2V2","UBE2N","H2AFX","CHAF1A","SETMAR","RECQL4","MPLKIP","DCLRE1A","DCLRE1B","PRPF19","RECQL","RECQL5","HELQ","RDM1","NABP2","ATRIP","MDC1","RAD1","RAD9A","HUS1","RAD17","TP53","TP53BP1","TOPBP1","CLK2","PER1"))

Extract DNA repair genes from GnomAD data by repair pathway

GnomAD_stats_LOF_VEP_ENSTcanonical_BER <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","XRCC1","PCNA","PARP1","PARP2","APEX1","FEN1","POLE","POLD1","LIG3","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","LIG1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_MMR <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","POLD1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_NER <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("RPA1","RPA2","RPA3","RPA4","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC1","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_HR <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","PALB2","MUS81","EME1","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_FA <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_DR <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("MGMT","ALKBH2","ALKBH3")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("XRCC1","PARP1","PARP2","APEX1","FEN1","OGG1","UNG","SMUG1","MBD4","TDG","MUTYH","NTHL1","MPG","NEIL1","NEIL2","NEIL3","APEX2","PNKP","APLF","PARP3","MSH3","MSH4","MSH5","MLH3","PMS1","PMS2P3","EXO1","RFC1","RFC2","RFC3","RFC4","RFC5","PCNA","ERCC4","ERCC2","ERCC5","XPC","ERCC6","GTF2H2","ERCC3","XPA","RAD23B","POLE","POLD1","RAD23A","LIG3","CETN2","DDB1","DDB2","GTF2H1","GTF2H3","GTF2H4","GTF2H5","CDK7","CCNH","MNAT1","LIG1","ERCC8","UVSSA","XAB2","MMS19","ATM","RPA1","RPA2","RPA3","RPA4","RAD51","NBN","RAD50","CHEK2","MRE11A","RAD52","RAD51B","DMC1","XRCC2","XRCC3","RAD54L","RAD54B","SHFM1","RBBP8","SLX1A","SLX1B","GEN1","ATR","ERCC1","FANCA","FANCB","FANCC","FANCD2","FANCE","FANCF","FANCG","FANCI","FANCL","FANCM","PALB2","CHEK1","SLX4","FAN1","C1orf86"," C19orf40","MUS81","EME1","XRCC6","XRCC5","PRKDC","LIG4","XRCC4","DCLRE1C","NHEJ1","MGMT","ALKBH2","ALKBH3")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER <- filter(GnomAD_stats_LOF_VEP,SYMBOL%in%c("PAXIP1","BLM","KMT2C","CRIP1","CDK12","BAP1","BARD1","WRN","BUB1","CENPE","ZW10","TTK","KNTC1","AURKB","POLB","POLH","POLQ","TDP1","TDP2","NUDT1","DUT","RRM2B","POLG","REV3L","MAD2L2","REV1","POLI","POLK","POLL","POLM","POLN","TREX1","TREX2","APTX","SPO11","ENDOV","UBE2A","UBE2B","RAD18","SHPRH","HLTF","RNF168","SPRTN","RNF8","RNF4","UBE2V2","UBE2N","H2AFX","CHAF1A","SETMAR","RECQL4","MPLKIP","DCLRE1A","DCLRE1B","PRPF19","RECQL","RECQL5","HELQ","RDM1","NABP2","ATRIP","MDC1","RAD1","RAD9A","HUS1","RAD17","TP53","TP53BP1","TOPBP1","CLK2","PER1")) %>% filter(CANONICAL=="YES")
GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair_OTHER <- bind_rows(
  GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair,
  GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER)

Calculate sample and GnomAD totals for DNA repair genes and extract into new data frame

a <- c("BER","MMR","NER","HR","FA","NHEJ","DR","ALL","OTHER")
b <- c(sum(genes_BER$Total_Gene_Count),
       sum(genes_MMR$Total_Gene_Count),
       sum(genes_NER$Total_Gene_Count),
       sum(genes_HR$Total_Gene_Count),
       sum(genes_FA$Total_Gene_Count),
       sum(genes_NHEJ$Total_Gene_Count),
       sum(genes_DR$Total_Gene_Count),
       sum(genes_ALLrepair$Total_Gene_Count),
       sum(genes_OTHER$Total_Gene_Count))
c <- c((n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_BER)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_MMR)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_NER)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_HR)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_FA)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_DR)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair)),
       (n_distinct(masterfile_eoc_goodQ2$Sample)*2*nrow(GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER)))
d <- c(sum(GnomAD_stats_LOF_VEP_ENSTcanonical_BER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_MMR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_HR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_FA$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_DR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair$FILTER_RF_LOF_HIGHIMPACT_AC_0.005),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005))
e <- c(sum(GnomAD_stats_LOF_VEP_ENSTcanonical_BER$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_MMR$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NER$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_HR$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_FA$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_DR$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair$MAX_AN),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER$MAX_AN))
f <- c(sum(GnomAD_stats_LOF_VEP_ENSTcanonical_BER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_MMR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_HR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_FA$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_DR$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE))
g <- c(sum(GnomAD_stats_LOF_VEP_ENSTcanonical_BER$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_MMR$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NER$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_HR$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_FA$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_NHEJ$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_DR$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_ALLrepair$MAX_AN_NFE),
       sum(GnomAD_stats_LOF_VEP_ENSTcanonical_OTHER$MAX_AN_NFE))

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary <- tibble("DNA_Repair_Pathway"=a,
                                                                                                            "Sample_Variant_Count"=b,
                                                                                                            "Total_Number_of_Sample_Alleles"=c,
                                                                                                            "GnomAD_Variant_Count"=d,
                                                                                                            "Total_Number_of_GnomAD_Variant_Alleles"=e,
                                                                                                            "GnomAD_Variant_Count_NFE"=f,
                                                                                                            "Total_Number_of_GnomAD_Variant_Alleles_NFE"=g)
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2 <-  add_row(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary,
        DNA_Repair_Pathway="TOTAL",
        Sample_Variant_Count=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,2]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,2]),
        Total_Number_of_Sample_Alleles=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,3]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,3]),
        GnomAD_Variant_Count=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,4]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,4]),
        Total_Number_of_GnomAD_Variant_Alleles=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,5]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,5]),
        GnomAD_Variant_Count_NFE=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,6]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,6]),
        Total_Number_of_GnomAD_Variant_Alleles_NFE=as.numeric(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[8,7]+masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary[9,7]))

Calculate Fisher’s test values for DNA repair genes

fisherresults_DNArepair <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2[,c(2,3,4,5)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$P_value_Fisher_0.005 = sapply(fisherresults_DNArepair,function(x) round(x$p.value,25))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$OR_0.005 = sapply(fisherresults_DNArepair,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$`95%CI_0.005` = sapply(fisherresults_DNArepair,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_DNArepair_NFE <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2[,c(2,3,6,7)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$P_value_Fisher_0.005_NFE = sapply(fisherresults_DNArepair_NFE,function(x) round(x$p.value,25))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$OR_0.005_NFE = sapply(fisherresults_DNArepair_NFE,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2$`95%CI_0.005_NFE` = sapply(fisherresults_DNArepair_NFE,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_DNArepair_other <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other[,c(221,229,243,244)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$P_value_Fisher_0.005 = sapply(fisherresults_DNArepair_other,function(x) round(x$p.value,12))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$OR_0.005 = sapply(fisherresults_DNArepair_other,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$`95%CI_0.005` = sapply(fisherresults_DNArepair_other,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

fisherresults_DNArepair_other_NFE <- apply(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other[,c(221,232,243,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$P_value_Fisher_0.005_NFE = sapply(fisherresults_DNArepair_other_NFE,function(x) round(x$p.value,12))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$OR_0.005_NFE = sapply(fisherresults_DNArepair_other_NFE,function(x) round(x$estimate,4))
masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other$`95%CI_0.005_NFE` = sapply(fisherresults_DNArepair_other_NFE,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))

Remove redundant data objects and output variants and summary results table for DNA repair genes

rm(list=(ls(pattern="^genes")))
rm(list=(ls(pattern="^variants")))
rm(list=(ls(pattern="^fisherresults")))

write_excel_csv(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other,path="masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_other.csv",na=".",append=FALSE,col_names=TRUE)
write_excel_csv(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_popmax0.005_DNArepair_summary2,path="masterfile_eoc_DNArepair_fisherresults.csv",na=".",append=FALSE,col_names=TRUE)