Thesis Missense Variant AF Analysis Script

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

Load tidyverse package

library("tidyverse")

Import raw vcf data from “masterfile” and exclude variants exceeding specified GnomAD non-cancer AF

masterfile <- 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)
  filter(GnomAD_v2.1_non_cancer_AF<=0.001)
  # filter(GnomAD_v2.1_non_cancer_AF<=0.0005)
  # filter(GnomAD_v2.1_non_cancer_AF<=0.0001)

Exclude non-epithelial, uterine-only, BRCA +ve and other mutation +ve samples

ViP_Complete_Cohort <- read.delim("ViP_LoF_Complete_Cohort.txt", stringsAsFactors=FALSE)
ViP_Complete_Cohort_list <- ViP_Complete_Cohort[,1]

masterfile_eoc <- filter(masterfile,Sample%in%c(ViP_Complete_Cohort_list))
rm(masterfile)

Append patient path data

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

Exclude low-grade tumour samples

masterfile_eoc2 <- filter(masterfile_eoc_withPath,(Histopath.Type!="Borderline tumour")&
                            (Histopath.Type!="Mucinous")&
                            (Histopath.Type!="Clear cell")&
                            (Histopath.Type!="Low-grade serous")&
                            (Histopath.Type!="Low-grade endometrioid")&
                            (Histopath.Type!="Serous (?low-grade/borderline)")&
                            (Histopath.Type!="Mixed Mullerian")&
                            (Histopath.Type!="Mixed EAOC"))
rm(masterfile_eoc_withPath)

Exclude low-quality variants, GnomAD RF-flagged variants and non-protein-coding genes

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)) %>% 
  filter(BIOTYPE%in%c("protein_coding"))
Total_Sample_Alleles <- (n_distinct(masterfile_eoc_goodQ$Sample)*2)
rm(masterfile_eoc2)

Exclude HIGH impact variants (plus variants with Consequence == “protein_altering_variant”)

masterfile_eoc_goodQ_biotype_MS <- filter(masterfile_eoc_goodQ, IMPACT=="MODERATE") %>% 
  filter(str_detect(Consequence,"missense_variant"))

Exclude NM canonical (RefSeq) variants and convert MPC strings to integers/doubles

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical <- filter(masterfile_eoc_goodQ_biotype_MS, CANONICAL=="YES") 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical$MPC <- as.numeric(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical$MPC) %>% replace_na(0)
rm(masterfile_eoc_goodQ_biotype_MS)

Divide data into different prediction score groups

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (MPC<2))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (MPC>=2))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (MPC>=2.5))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (MPC>=3))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED<15))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=15))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=20))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=25))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=30))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score<0.45))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.45))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.5))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.55))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.6))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.7))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (REVEL_score>=0.8))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=15)&(REVEL_score>=0.45))  
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical, (CADD_PHRED>=15)&(REVEL_score>=0.45)&(MPC>=2))
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8 <- bind_rows(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8,.id=NULL) %>% 
  distinct(.keep_all=TRUE)

rm(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical)

Sample variant counts and AFs

variant_counts_AFs <- function(vcf){
variants_heterozygous <- vcf %>% 
  select(HGVSc,Sample.GT) %>% 
  group_by(HGVSc,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT%in%c("'0/1","'1/0"))) %>% 
  dplyr::rename(alleles = "n()") %>% 
  group_by(HGVSc) %>% 
  summarise(sum(alleles))
variants_homozygous <- vcf %>% 
  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/Total_Sample_Alleles)
vcf <- left_join(vcf, variants[,c(1,4:5)],by="HGVSc",copy=FALSE)
}

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2  <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8 <- variant_counts_AFs(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8)

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

Output list of genes and variants with sample AF > 0.01

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_genesandvariants0.01 <- bind_rows(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8) %>% 
  filter(Sample_AF > 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_genesandvariants0.01[,c(2:6,25:54)] %>% 
  distinct(.keep_all = FALSE) %>% 
  write_excel_csv(path="masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_genesandvariants0.01.csv",na=".",append=FALSE,col_names=TRUE)
rm(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_genesandvariants0.01)

Exclude variants with sample AF > 0.01̛

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3,Sample_AF < 0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30,Sample_AF < 0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8,Sample_AF < 0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45,Sample_AF < 0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01 <- filter(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8,Sample_AF < 0.01)

rm(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45,masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8)

Sample gene counts and frequencies

gene_counts_AF <- function(vcf){
genes_oneVarAllele <- vcf %>% select(SYMBOL,Sample.GT) %>% 
  select(SYMBOL,Sample.GT) %>% 
  group_by(SYMBOL,Sample.GT) %>% 
  summarise(n()) %>% 
  filter((Sample.GT%in%c("'0/1","'1/0"))) %>% 
  dplyr::rename(gene_Var = "n()") %>% 
  group_by(SYMBOL) %>% 
  summarise(sum(gene_Var))
genes_twoVarAllele <- vcf %>% select(SYMBOL,Sample.GT) %>% 
  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/Total_Sample_Alleles)
vcf <- left_join(vcf, genes[,c(1,4:5)],by="SYMBOL",copy=FALSE)
}

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01 <- gene_counts_AF(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01)

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

Add GnomAD gene-level data

GnomAD_stats_MS <- read.delim("GnomAD.gene.stats.v2.1_missense_AF_0.001.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  dplyr::rename(SYMBOL=X.Symbol) 
GnomAD_stats_MS <- GnomAD_stats_MS[,c(1:13,156:281)] 

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,14:19)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,20:25)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,26:31)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,32:37)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,38:43)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,44:49)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,50:55)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,56:61)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,62:67)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,68:73)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,74:79)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,80:85)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,86:91)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,92:97)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,104:109)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,116:121)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,134:139)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,122:127)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_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)))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01,GnomAD_stats_MS[,c(2,6:13,128:133)],by="Feature",copy=FALSE) 
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01_withGnomADstats <-
mutate_at(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_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)))

rm(list=(ls(pattern="_sampleAF0.01$")))

Create data frames with Agilent SureSelect whole exome genes (ENST and protein-coding) and MPC-only transcripts

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_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"))
MPC_transcripts <- read.delim("MPC_transcripts_with_canonicity.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  filter(CANONICAL%in%c("YES")) %>% 
  select(Feature,CANONICAL)

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, "AgilentSSv6"="YES") %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)

allGenes_AgilentSSv6_list_protein_coding_MPConly<- 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, "AgilentSSv6"="YES") %>% 
  distinct(.keep_all=TRUE) %>% 
  left_join(MPC_transcripts,by="Feature",copy=FALSE) %>% 
  filter(CANONICAL%in%c("YES"))

allGenes_AgilentSSv6_list_protein_coding_only2 <- select(allGenes_AgilentSSv6_list_protein_coding_only,Feature,AgilentSSv6)

Sum ACs across sample/GnomAD for each prediction score class (MPC figures include transcripts with MPC scores only)

GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding <- filter(GnomAD_stats_MS,CANONICAL%in%c("YES")) %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  left_join(allGenes_AgilentSSv6_list_protein_coding_only2,by="Feature",copy=FALSE) %>% 
  filter(AgilentSSv6%in%c("YES")) %>% 
  distinct(.keep_all=TRUE)
GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC <- filter(GnomAD_stats_MS,CANONICAL%in%c("YES")) %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  left_join(allGenes_AgilentSSv6_list_protein_coding_only2,by="Feature",copy=FALSE) %>% 
  filter(AgilentSSv6%in%c("YES")) %>% 
  left_join(MPC_transcripts,by="Feature",copy=FALSE) %>% 
  filter(CANONICAL.y%in%c("YES")) %>% 
  distinct(.keep_all=TRUE)

GnomAD_stats_MS_MAX_AN_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,6]) %>% as.numeric()
GnomAD_stats_MS_MAX_AN_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,10]) %>% as.numeric()
GnomAD_stats_MS_MAX_AN_MPC_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,6]) %>% as.numeric()
GnomAD_stats_MS_MAX_AN_MPC_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,10]) %>% as.numeric()
GnomAD_stats_MS_MPC.2_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,14]) %>% as.numeric()
GnomAD_stats_MS_MPC.2_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,17]) %>% as.numeric()
GnomAD_stats_MS_MPC2_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,20]) %>% as.numeric()
GnomAD_stats_MS_MPC2_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,23]) %>% as.numeric()
GnomAD_stats_MS_MPC2.5_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,26]) %>% as.numeric()
GnomAD_stats_MS_MPC2.5_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,29]) %>% as.numeric()
GnomAD_stats_MS_MPC3_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,32]) %>% as.numeric()
GnomAD_stats_MS_MPC3_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,35]) %>% as.numeric()
GnomAD_stats_MS_CADD.15_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,38]) %>% as.numeric()
GnomAD_stats_MS_CADD.15_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,41]) %>% as.numeric()
GnomAD_stats_MS_CADD15_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,44]) %>% as.numeric()
GnomAD_stats_MS_CADD15_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,47]) %>% as.numeric()
GnomAD_stats_MS_CADD20_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,50]) %>% as.numeric()
GnomAD_stats_MS_CADD20_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,53]) %>% as.numeric()
GnomAD_stats_MS_CADD25_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,56]) %>% as.numeric()
GnomAD_stats_MS_CADD25_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,59]) %>% as.numeric()
GnomAD_stats_MS_CADD30_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,62]) %>% as.numeric()
GnomAD_stats_MS_CADD30_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,65]) %>% as.numeric()
GnomAD_stats_MS_REVEL.0.45_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,68]) %>% as.numeric()
GnomAD_stats_MS_REVEL.0.45_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,71]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.45_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,74]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.45_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,77]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.5_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,80]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.5_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,83]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.55_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,86]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.55_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,89]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.6_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,92]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.6_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,95]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.7_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,104]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.7_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,107]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.8_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,116]) %>% as.numeric()
GnomAD_stats_MS_REVEL0.8_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,119]) %>% as.numeric()
GnomAD_stats_MS_CADD15REVEL0.45_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,134]) %>% as.numeric()
GnomAD_stats_MS_CADD15REVEL0.45_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding[,137]) %>% as.numeric()
GnomAD_stats_MS_MPC2CADD15REVEL0.45_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,122]) %>% as.numeric()
GnomAD_stats_MS_MPC2CADD15REVEL0.45_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,125]) %>% as.numeric()
GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,128]) %>% as.numeric()
GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE_total <- sum(GnomAD_stats_MS_ENSTcanonical_AgilentSSv6_protein_coding_MPC[,131]) %>% as.numeric()

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total <- Total_Sample_Alleles*(nrow(allGenes_AgilentSSv6_list_protein_coding_only)) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total <- Total_Sample_Alleles*(nrow(allGenes_AgilentSSv6_list_protein_coding_MPConly)) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total <- masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_sampleAF0.01_withGnomADstats[,c(31,36,37,221)] %>% distinct(.keep_all=TRUE) %>% select(4) %>% colSums(na.rm=TRUE) %>% as.numeric()

Calculate odds ratios using matrices

MPC.2 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_total,
                 masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_total,
                 GnomAD_stats_MS_MPC.2_total,
                 GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC.2_total),
               nrow=2,byrow=TRUE)
MPC.2_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_total,
                     masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC.2_total,
                     GnomAD_stats_MS_MPC.2_NFE_total,
                     GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC.2_NFE_total),
                   nrow=2,byrow=TRUE)
MPC2 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_total,
                     masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_total,
                     GnomAD_stats_MS_MPC2_total,
                     GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC2_total),
                   nrow=2,byrow=TRUE)
MPC2_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_total,
                         masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_total,
                         GnomAD_stats_MS_MPC2_NFE_total,
                         GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC2_NFE_total),
                       nrow=2,byrow=TRUE)
MPC2.5 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_total,
                     masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_total,
                     GnomAD_stats_MS_MPC2.5_total,
                     GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC2.5_total),
                   nrow=2,byrow=TRUE)
MPC2.5_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_total,
                         masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2.5_total,
                         GnomAD_stats_MS_MPC2.5_NFE_total,
                         GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC2.5_NFE_total),
                       nrow=2,byrow=TRUE)
MPC3 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_total,
                 masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_total,
                 GnomAD_stats_MS_MPC3_total,
                 GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC3_total),
               nrow=2,byrow=TRUE)
MPC3_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_total,
                     masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC3_total,
                     GnomAD_stats_MS_MPC3_NFE_total,
                     GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC3_NFE_total),
                   nrow=2,byrow=TRUE)
CADD.15 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_total,
                   masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_total,
                   GnomAD_stats_MS_CADD.15_total,
                   GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD.15_total),
                 nrow=2,byrow=TRUE)
CADD.15_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_total,
                       masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD.15_total,
                       GnomAD_stats_MS_CADD.15_NFE_total,
                       GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD.15_NFE_total),
                     nrow=2,byrow=TRUE)
CADD15 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_total,
                      masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_total,
                      GnomAD_stats_MS_CADD15_total,
                      GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD15_total),
                    nrow=2,byrow=TRUE)
CADD15_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15_total,
                          GnomAD_stats_MS_CADD15_NFE_total,
                          GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD15_NFE_total),
                        nrow=2,byrow=TRUE)
CADD20 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_total,
                      masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_total,
                      GnomAD_stats_MS_CADD20_total,
                      GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD20_total),
                    nrow=2,byrow=TRUE)
CADD20_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD20_total,
                          GnomAD_stats_MS_CADD20_NFE_total,
                          GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD20_NFE_total),
                        nrow=2,byrow=TRUE)
CADD25 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_total,
                      masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_total,
                      GnomAD_stats_MS_CADD25_total,
                      GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD25_total),
                    nrow=2,byrow=TRUE)
CADD25_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD25_total,
                          GnomAD_stats_MS_CADD25_NFE_total,
                          GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD25_NFE_total),
                        nrow=2,byrow=TRUE)
CADD30 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_total,
                   masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_total,
                   GnomAD_stats_MS_CADD30_total,
                   GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD30_total),
                 nrow=2,byrow=TRUE)
CADD30_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_total,
                       masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD30_total,
                       GnomAD_stats_MS_CADD30_NFE_total,
                       GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD30_NFE_total),
                     nrow=2,byrow=TRUE)
REVEL.0.45 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_total,
                      masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_total,
                      GnomAD_stats_MS_REVEL.0.45_total,
                      GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL.0.45_total),
                    nrow=2,byrow=TRUE)
REVEL.0.45_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL.0.45_total,
                          GnomAD_stats_MS_REVEL.0.45_NFE_total,
                          GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL.0.45_NFE_total),
                        nrow=2,byrow=TRUE)
REVEL0.45 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_total,
                          GnomAD_stats_MS_REVEL0.45_total,
                          GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.45_total),
                        nrow=2,byrow=TRUE)
REVEL0.45_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_total,
                              masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_total,
                              GnomAD_stats_MS_REVEL0.45_NFE_total,
                              GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.45_NFE_total),
                            nrow=2,byrow=TRUE)
REVEL0.5 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_total,
                          GnomAD_stats_MS_REVEL0.5_total,
                          GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.5_total),
                        nrow=2,byrow=TRUE)
REVEL0.5_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_total,
                              masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.5_total,
                              GnomAD_stats_MS_REVEL0.5_NFE_total,
                              GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.5_NFE_total),
                            nrow=2,byrow=TRUE)
REVEL0.55 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_total,
                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_total,
                          GnomAD_stats_MS_REVEL0.55_total,
                          GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.55_total),
                        nrow=2,byrow=TRUE)
REVEL0.55_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_total,
                              masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.55_total,
                              GnomAD_stats_MS_REVEL0.55_NFE_total,
                              GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.55_NFE_total),nrow=2,byrow=TRUE)
REVEL0.6 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_total,
                         masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_total,
                         GnomAD_stats_MS_REVEL0.6_total,
                         GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.6_total),
                       nrow=2,byrow=TRUE)
REVEL0.6_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_total,
                             masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.6_total,
                             GnomAD_stats_MS_REVEL0.6_NFE_total,
                             GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.6_NFE_total),
                           nrow=2,byrow=TRUE)
REVEL0.7 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_total,
                         masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_total,
                         GnomAD_stats_MS_REVEL0.7_total,
                         GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.7_total),
                       nrow=2,byrow=TRUE)
REVEL0.7_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_total,
                             masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.7_total,
                             GnomAD_stats_MS_REVEL0.7_NFE_total,
                             GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.7_NFE_total),
                           nrow=2,byrow=TRUE)
REVEL0.8 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_total,
                     masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_total,
                     GnomAD_stats_MS_REVEL0.8_total,
                     GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_REVEL0.8_total),
                   nrow=2,byrow=TRUE)
REVEL0.8_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_total,
                         masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.8_total,
                         GnomAD_stats_MS_REVEL0.8_NFE_total,
                         GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_REVEL0.8_NFE_total),
                       nrow=2,byrow=TRUE)
CADD15REVEL0.45 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_total,
                            masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_total,
                            GnomAD_stats_MS_CADD15REVEL0.45_total,
                            GnomAD_stats_MS_MAX_AN_total-GnomAD_stats_MS_CADD15REVEL0.45_total),
                          nrow=2,byrow=TRUE)
CADD15REVEL0.45_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_total,
                                masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_CADD15REVEL0.45_total,
                                GnomAD_stats_MS_CADD15REVEL0.45_NFE_total,
                                GnomAD_stats_MS_MAX_AN_NFE_total-GnomAD_stats_MS_CADD15REVEL0.45_NFE_total),
                              nrow=2,byrow=TRUE)
MPC2CADD15REVEL0.45 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_total,
                                masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_total,
                                GnomAD_stats_MS_MPC2CADD15REVEL0.45_total,
                                GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC2CADD15REVEL0.45_total),
                              nrow=2,byrow=TRUE)
MPC2CADD15REVEL0.45_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_total,
                                    masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2CADD15REVEL0.45_total,
                                    GnomAD_stats_MS_MPC2CADD15REVEL0.45_NFE_total,
                                    GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC2CADD15REVEL0.45_NFE_total),nrow=2,byrow=TRUE)
MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8 <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total,
                                                          masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total,
                                                          GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total,
                                                          GnomAD_stats_MS_MAX_AN_MPC_total-GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total),
                                                        nrow=2,byrow=TRUE)
MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE <- matrix(c(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total,
                                                              masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_all_alleles_MPConly_total-masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_total,
                                                              GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE_total,
                                                              GnomAD_stats_MS_MAX_AN_MPC_NFE_total-GnomAD_stats_MS_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE_total),
                                                            nrow=2,byrow=TRUE)

oddsratio <- function (a, b = NULL, c = NULL, d = NULL, conf.level = 0.95, 
    p.calc.by.independence = TRUE) 
{
    if (is.matrix(a)) {
        if ((dim(a)[1] != 2L) | (dim(a)[2] != 2L)) {
            stop("Input matrix must be a 2x2 table.")
        }
        .a <- a[1, 1]
        .b <- a[1, 2]
        .c <- a[2, 1]
        .d <- a[2, 2]
        .data.name <- deparse(substitute(a))
    }
    else {
        .a <- a
        .b <- b
        .c <- c
        .d <- d
        .data.name <- paste(deparse(substitute(a)), deparse(substitute(b)), 
            deparse(substitute(c)), deparse(substitute(d)))
    }
    .MAT <- matrix(c(.a, .b, M1 <- .a + .b, .c, .d, M0 <- .c + 
        .d, N1 <- .a + .c, N0 <- .b + .d, Total <- .a + .b + 
        .c + .d), 3, 3)
    colnames(.MAT) <- c("Sample", "GnomAD", "Total") #("Disease", "Nondisease", "Total")
    rownames(.MAT) <- c("MS", "No MS", "Total") #("Exposed", "Nonexposed", "Total")
    class(.MAT) <- "table"
    print(.MAT)
    ESTIMATE <- (.a /.b)/(.c/.d)
    norm.pp <- qnorm(1 - (1 - conf.level)/2)
    if (p.calc.by.independence) {
        p.v <- 2 * (1 - pnorm(abs((.a - N1 * M1/Total)/sqrt(N1 * 
            N0 * M1 * M0/Total/Total/(Total - 1)))))
    }
    else {
        p.v <- 2 * (1 - pnorm(log(ifelse(ESTIMATE > 1, ESTIMATE, 
            1/ESTIMATE))/sqrt(1/.a + 1/.b + 1/.c + 1/.d)))
    }
    ORL <- ESTIMATE * exp(-norm.pp * sqrt(1/.a + 1/.b + 1/.c + 
        1/.d))
    ORU <- ESTIMATE * exp(norm.pp * sqrt(1/.a + 1/.b + 1/.c + 
        1/.d)) %>% signif(digits=7)
    CINT <- paste(signif(ORL,digits = 7),signif(ORU,digits = 7),sep="~")
    attr(CINT, "conf.level") <- conf.level
    RVAL <- list(p.value = p.v, conf.int = CINT, estimate = ESTIMATE, 
        method = "Odds ratio estimate and its significance probability", 
        data.name = .data.name)
    class(RVAL) <- "htest"
    return(RVAL)
}

OR_MPC.2 <- oddsratio(MPC.2,p.calc.by.independence=TRUE)
OR_MPC.2_NFE <- oddsratio(MPC.2_NFE,p.calc.by.independence=TRUE)
OR_MPC2 <- oddsratio(MPC2,p.calc.by.independence=TRUE)
OR_MPC2_NFE <- oddsratio(MPC2_NFE,p.calc.by.independence=TRUE)
OR_MPC2.5 <- oddsratio(MPC2.5,p.calc.by.independence=TRUE)
OR_MPC2.5_NFE <- oddsratio(MPC2.5_NFE,p.calc.by.independence=TRUE)
OR_MPC3 <- oddsratio(MPC3,p.calc.by.independence=TRUE)
OR_MPC3_NFE <- oddsratio(MPC3_NFE,p.calc.by.independence=TRUE)
OR_CADD.15 <- oddsratio(CADD.15,p.calc.by.independence=TRUE)
OR_CADD.15_NFE <- oddsratio(CADD.15_NFE,p.calc.by.independence=TRUE)
OR_CADD15 <- oddsratio(CADD15,p.calc.by.independence=TRUE)
OR_CADD15_NFE <- oddsratio(CADD15_NFE,p.calc.by.independence=TRUE)
OR_CADD20 <- oddsratio(CADD20,p.calc.by.independence=TRUE)
OR_CADD20_NFE <- oddsratio(CADD20_NFE,p.calc.by.independence=TRUE)
OR_CADD25 <- oddsratio(CADD25,p.calc.by.independence=TRUE)
OR_CADD25_NFE <- oddsratio(CADD25_NFE,p.calc.by.independence=TRUE)
OR_CADD30 <- oddsratio(CADD30,p.calc.by.independence=TRUE)
OR_CADD30_NFE <- oddsratio(CADD30_NFE,p.calc.by.independence=TRUE)
OR_REVEL.0.45 <- oddsratio(REVEL.0.45,p.calc.by.independence=TRUE)
OR_REVEL.0.45_NFE <- oddsratio(REVEL.0.45_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.45 <- oddsratio(REVEL0.45,p.calc.by.independence=TRUE)
OR_REVEL0.45_NFE <- oddsratio(REVEL0.45_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.5 <- oddsratio(REVEL0.5,p.calc.by.independence=TRUE)
OR_REVEL0.5_NFE <- oddsratio(REVEL0.5_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.55 <- oddsratio(REVEL0.55,p.calc.by.independence=TRUE)
OR_REVEL0.55_NFE <- oddsratio(REVEL0.55_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.6 <- oddsratio(REVEL0.6,p.calc.by.independence=TRUE)
OR_REVEL0.6_NFE <- oddsratio(REVEL0.6_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.7 <- oddsratio(REVEL0.7,p.calc.by.independence=TRUE)
OR_REVEL0.7_NFE <- oddsratio(REVEL0.7_NFE,p.calc.by.independence=TRUE)
OR_REVEL0.8 <- oddsratio(REVEL0.8,p.calc.by.independence=TRUE)
OR_REVEL0.8_NFE <- oddsratio(REVEL0.8_NFE,p.calc.by.independence=TRUE)
OR_CADD15REVEL0.45 <- oddsratio(CADD15REVEL0.45,p.calc.by.independence=TRUE)
OR_CADD15REVEL0.45_NFE <- oddsratio(CADD15REVEL0.45_NFE,p.calc.by.independence=TRUE)
OR_MPC2CADD15REVEL0.45 <- oddsratio(MPC2CADD15REVEL0.45,p.calc.by.independence=TRUE)
OR_MPC2CADD15REVEL0.45_NFE <- oddsratio(MPC2CADD15REVEL0.45_NFE,p.calc.by.independence=TRUE)
OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8 <- oddsratio(MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8,p.calc.by.independence=TRUE)
OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE <- oddsratio(MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE,p.calc.by.independence=TRUE)

Combine all results into one spreadsheet, and save

a = c(OR_MPC.2$data.name,OR_MPC.2_NFE$data.name,OR_MPC2$data.name,OR_MPC2_NFE$data.name,OR_MPC2.5$data.name,OR_MPC2.5_NFE$data.name,OR_MPC3$data.name,OR_MPC3_NFE$data.name,OR_CADD.15$data.name,OR_CADD.15_NFE$data.name,OR_CADD15$data.name,OR_CADD15_NFE$data.name,OR_CADD20$data.name,OR_CADD20_NFE$data.name,OR_CADD25$data.name,OR_CADD25_NFE$data.name,OR_CADD30$data.name,OR_CADD30_NFE$data.name,OR_REVEL.0.45$data.name,OR_REVEL.0.45_NFE$data.name,OR_REVEL0.45$data.name,OR_REVEL0.45_NFE$data.name,OR_REVEL0.5$data.name,OR_REVEL0.5_NFE$data.name,OR_REVEL0.55$data.name,OR_REVEL0.55_NFE$data.name,OR_REVEL0.6$data.name,OR_REVEL0.6_NFE$data.name,OR_REVEL0.7$data.name,OR_REVEL0.7_NFE$data.name,OR_REVEL0.8$data.name,OR_REVEL0.8_NFE$data.name,OR_CADD15REVEL0.45$data.name,OR_CADD15REVEL0.45_NFE$data.name,OR_MPC2CADD15REVEL0.45$data.name,OR_MPC2CADD15REVEL0.45_NFE$data.name,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8$data.name,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE$data.name)
b = c(OR_MPC.2$estimate,OR_MPC.2_NFE$estimate,OR_MPC2$estimate,OR_MPC2_NFE$estimate,OR_MPC2.5$estimate,OR_MPC2.5_NFE$estimate,OR_MPC3$estimate,OR_MPC3_NFE$estimate,OR_CADD.15$estimate,OR_CADD.15_NFE$estimate,OR_CADD15$estimate,OR_CADD15_NFE$estimate,OR_CADD20$estimate,OR_CADD20_NFE$estimate,OR_CADD25$estimate,OR_CADD25_NFE$estimate,OR_CADD30$estimate,OR_CADD30_NFE$estimate,OR_REVEL.0.45$estimate,OR_REVEL.0.45_NFE$estimate,OR_REVEL0.45$estimate,OR_REVEL0.45_NFE$estimate,OR_REVEL0.5$estimate,OR_REVEL0.5_NFE$estimate,OR_REVEL0.55$estimate,OR_REVEL0.55_NFE$estimate,OR_REVEL0.6$estimate,OR_REVEL0.6_NFE$estimate,OR_REVEL0.7$estimate,OR_REVEL0.7_NFE$estimate,OR_REVEL0.8$estimate,OR_REVEL0.8_NFE$estimate,OR_CADD15REVEL0.45$estimate,OR_CADD15REVEL0.45_NFE$estimate,OR_MPC2CADD15REVEL0.45$estimate,OR_MPC2CADD15REVEL0.45_NFE$estimate,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8$estimate,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE$estimate)
c = c(OR_MPC.2$conf.int,OR_MPC.2_NFE$conf.int,OR_MPC2$conf.int,OR_MPC2_NFE$conf.int,OR_MPC2.5$conf.int,OR_MPC2.5_NFE$conf.int,OR_MPC3$conf.int,OR_MPC3_NFE$conf.int,OR_CADD.15$conf.int,OR_CADD.15_NFE$conf.int,OR_CADD15$conf.int,OR_CADD15_NFE$conf.int,OR_CADD20$conf.int,OR_CADD20_NFE$conf.int,OR_CADD25$conf.int,OR_CADD25_NFE$conf.int,OR_CADD30$conf.int,OR_CADD30_NFE$conf.int,OR_REVEL.0.45$conf.int,OR_REVEL.0.45_NFE$conf.int,OR_REVEL0.45$conf.int,OR_REVEL0.45_NFE$conf.int,OR_REVEL0.5$conf.int,OR_REVEL0.5_NFE$conf.int,OR_REVEL0.55$conf.int,OR_REVEL0.55_NFE$conf.int,OR_REVEL0.6$conf.int,OR_REVEL0.6_NFE$conf.int,OR_REVEL0.7$conf.int,OR_REVEL0.7_NFE$conf.int,OR_REVEL0.8$conf.int,OR_REVEL0.8_NFE$conf.int,OR_CADD15REVEL0.45$conf.int,OR_CADD15REVEL0.45_NFE$conf.int,OR_MPC2CADD15REVEL0.45$conf.int,OR_MPC2CADD15REVEL0.45_NFE$conf.int,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8$conf.int,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE$conf.int)
d = c(OR_MPC.2$p.value,OR_MPC.2_NFE$p.value,OR_MPC2$p.value,OR_MPC2_NFE$p.value,OR_MPC2.5$p.value,OR_MPC2.5_NFE$p.value,OR_MPC3$p.value,OR_MPC3_NFE$p.value,OR_CADD.15$p.value,OR_CADD.15_NFE$p.value,OR_CADD15$p.value,OR_CADD15_NFE$p.value,OR_CADD20$p.value,OR_CADD20_NFE$p.value,OR_CADD25$p.value,OR_CADD25_NFE$p.value,OR_CADD30$p.value,OR_CADD30_NFE$p.value,OR_REVEL.0.45$p.value,OR_REVEL.0.45_NFE$p.value,OR_REVEL0.45$p.value,OR_REVEL0.45_NFE$p.value,OR_REVEL0.5$p.value,OR_REVEL0.5_NFE$p.value,OR_REVEL0.55$p.value,OR_REVEL0.55_NFE$p.value,OR_REVEL0.6$p.value,OR_REVEL0.6_NFE$p.value,OR_REVEL0.7$p.value,OR_REVEL0.7_NFE$p.value,OR_REVEL0.8$p.value,OR_REVEL0.8_NFE$p.value,OR_CADD15REVEL0.45$p.value,OR_CADD15REVEL0.45_NFE$p.value,OR_MPC2CADD15REVEL0.45$p.value,OR_MPC2CADD15REVEL0.45_NFE$p.value,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8$p.value,OR_MPC2_OR_CADD15REVEL0.45_OR_CADD30_OR_REVEL0.8_NFE$p.value)

ORs_masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_sample0.01 <- tibble(
  "Filtering"=a,
  "Odds_ratio"=b,
  "95%_CI"=c,
  "p_value"=d
)

write_tsv(ORs_masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_sample0.01, path="ORs_masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_0.001_sample0.01.tsv")

Calculate ratios and Fisher’s test values for genes containing enriched variants with corresponding prediction scores, and save output

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios <- 
  mutate(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats,
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_0.001"=Sample_Gene_Freq/FILTER_RF_MS_AF_0.001_REVEL0.45_ADJ,
         "Sample_Gene_LOF_Freq_Ratio_GnomAD_NFE_0.001"=Sample_Gene_Freq/FILTER_RF_MS_AF_0.001_REVEL0.45_NFE_ADJ)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios2 <- mutate(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_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)

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios2,Total_Case_Alleles=Total_Sample_Alleles)
masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_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_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_ENST <- apply(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(221,231,241,223)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

fisherresults2_ENST <- apply(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(221,234,241,227)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_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))

masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_REVEL0.45_sampleAF0.01_withGnomADstats_ratios_fisherresults3 %>% 
  write_excel_csv(path="masterfile_eoc_goodQ_biotype_MS_ENSTcanonical_0.001_REVEL0.45_sample0.01.csv",na=".",append=FALSE,col_names=TRUE)