Thesis LoF Variants Filtering and Analysis by Family History 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 with GnomAD non-cancer AF > 0.005

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)

Exclude non-epithelial, uterine-only and BRCA +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))

Append patient path data and family history info

ViP_OvCa_Path_Data <- read.delim("ViP_OvCa_Path_Data.txt", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  dplyr::rename("Sample"=Exome.ID)
NoFHxCa <- read.delim("NoFHxCa.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  select(Exome.ID,FHx) %>% 
  dplyr::rename("Sample"=Exome.ID)
masterfile_eoc_withPath_FHx <- left_join(masterfile_eoc,ViP_OvCa_Path_Data,by="Sample",copy=FALSE) %>% 
  left_join(NoFHxCa,by="Sample",copy=FALSE) %>% 
  mutate_at(vars("FHx"),funs(replace(., is.na(.), "YES")))
rm(ViP_OvCa_Path_Data,NoFHxCa)

Exclude low-grade tumour samples

masterfile_eoc2 <- filter(masterfile_eoc_withPath_FHx,(!Histopath.Type%in%c("Borderline tumour",
                                                                            "Mucinous",
                                                                            "Clear cell",
                                                                            "Low-grade serous",
                                                                            "Low-grade endometrioid",
                                                                            "Serous (?low-grade/borderline)",
                                                                            "Mixed Mullerian",
                                                                            "Mixed EAOC")))

Exclude samples with/without family history of BrOvCa in 1st and/or 2nd degree relatives

masterfile_eoc_FHx <- filter(masterfile_eoc2,FHx%in%c("YES"))

masterfile_eoc_noFHx <- filter(masterfile_eoc2,!FHx%in%c("YES"))

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

masterfile_eoc_FHx_goodQ <- filter(masterfile_eoc_FHx, (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))

masterfile_eoc_noFHx_goodQ <- filter(masterfile_eoc_noFHx, (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)

ViP_Discovery_Cohort <- read.delim("ViP_LoF_Discovery_Cohort.txt", stringsAsFactors=FALSE)
ViP_Discovery_Cohort_list <- ViP_Discovery_Cohort[,1]

masterfile_eoc_FHx_goodQ2 <- filter(masterfile_eoc_FHx_goodQ,Sample%in%c(ViP_Discovery_Cohort_list))
Total_Sample_Alleles_FHx <- (n_distinct(masterfile_eoc_FHx_goodQ2$Sample)*2)

masterfile_eoc_noFHx_goodQ2 <- filter(masterfile_eoc_noFHx_goodQ,Sample%in%c(ViP_Discovery_Cohort_list))
Total_Sample_Alleles_noFHx <- (n_distinct(masterfile_eoc_noFHx_goodQ2$Sample)*2)

Exclude MODERATE impact variants

masterfile_eoc_FHx_goodQ_HIGH <- filter(masterfile_eoc_FHx_goodQ, IMPACT=="HIGH")

masterfile_eoc_FHx_goodQ2_HIGH <- filter(masterfile_eoc_FHx_goodQ2, IMPACT=="HIGH")

masterfile_eoc_noFHx_goodQ_HIGH <- filter(masterfile_eoc_noFHx_goodQ, IMPACT=="HIGH")

masterfile_eoc_noFHx_goodQ2_HIGH <- filter(masterfile_eoc_noFHx_goodQ2, IMPACT=="HIGH")

Separate Ensembl canonical and RefSeq canonical variants

masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical <- filter(masterfile_eoc_FHx_goodQ_HIGH,CANONICAL=="YES") 

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical <- filter(masterfile_eoc_FHx_goodQ2_HIGH,CANONICAL=="YES") 

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical <- filter(masterfile_eoc_FHx_goodQ2_HIGH,Ensembl94_refseq_mrna!="is.na") %>% 
  filter(!Ensembl94_refseq_mrna_STATUS%in%c(".|.",
                                            ".|.|.",
                                            "MODEL",
                                            "INFERRED",
                                            "PREDICTED",
                                            "PROVISIONAL",
                                            "WGS"))

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical <- filter(masterfile_eoc_noFHx_goodQ_HIGH,CANONICAL=="YES") 

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical <- filter(masterfile_eoc_noFHx_goodQ2_HIGH,CANONICAL=="YES") 

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical <- filter(masterfile_eoc_noFHx_goodQ2_HIGH,Ensembl94_refseq_mrna!="is.na") %>% 
  filter(!Ensembl94_refseq_mrna_STATUS%in%c(".|.",
                                            ".|.|.",
                                            "MODEL",
                                            "INFERRED",
                                            "PREDICTED",
                                            "PROVISIONAL",
                                            "WGS"))

Sample variant counts and AFs

variant_counts_AFs <- function(vcf,Total_Sample_Alleles){
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_FHx_goodQ_HIGH_ENSTcanonical_2 <- variant_counts_AFs(masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical,Total_Sample_Alleles_FHx)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_2 <- variant_counts_AFs(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical,Total_Sample_Alleles_FHx)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_2 <- variant_counts_AFs(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical,Total_Sample_Alleles_FHx)

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_2 <- variant_counts_AFs(masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical,Total_Sample_Alleles_noFHx)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_2 <- variant_counts_AFs(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical,Total_Sample_Alleles_noFHx)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_2 <- variant_counts_AFs(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical,Total_Sample_Alleles_noFHx)

Exclude variants with sample AF >0.01̛

masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical_2,Sample_AF < 0.01)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_2,Sample_AF < 0.01)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01 <- filter(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_2,Sample_AF < 0.01)

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_2,Sample_AF < 0.01)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_2,Sample_AF < 0.01)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01 <- filter(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_2,Sample_AF < 0.01)

Output list of genes and variants with sample AF >0.01

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01 <- filter(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_2,Sample_AF > 0.01)
masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01[,c(2:6,25:54)] %>% 
  distinct(.keep_all = FALSE) %>% 
  write_excel_csv(path="masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01.csv",na=".",append=FALSE,col_names=TRUE)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01 <- filter(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_2,Sample_AF > 0.01)
masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01[,c(2:6,25:54)] %>% 
  distinct(.keep_all = FALSE) %>% 
  write_excel_csv(path="masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01.csv",na=".",append=FALSE,col_names=TRUE)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01 <- filter(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_2,Sample_AF > 0.01)
masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01[,c(2:6,25:54)] %>% 
  distinct(.keep_all = FALSE) %>% 
  write_excel_csv(path="masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_genesandvariants0.01.csv",na=".",append=FALSE,col_names=TRUE)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01 <- filter(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_2,Sample_AF > 0.01)
masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01[,c(2:6,25:54)] %>% 
  distinct(.keep_all = FALSE) %>% 
  write_excel_csv(path="masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_genesandvariants0.01.csv",na=".",append=FALSE,col_names=TRUE)

Sample gene counts and frequencies

gene_counts_AFs <- function(vcf,Total_Sample_Alleles){
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_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01,Total_Sample_Alleles_FHx)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01,Total_Sample_Alleles_FHx)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01,Total_Sample_Alleles_FHx)

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01,Total_Sample_Alleles_noFHx)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01,Total_Sample_Alleles_noFHx)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_2 <- gene_counts_AFs(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01,Total_Sample_Alleles_noFHx)

Add GnomAD gene-level data

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

GnomADD <- function(vcf,GnomAD){
  vcf <- left_join(vcf, GnomAD[,c(2,30:41)],by="Feature",copy=FALSE) 
  vcf <- mutate_at(vcf,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_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats <- GnomADD(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_2,GnomAD_stats_LOF_VEP)

Calculate Sample/GnomAD gene freq ratios for total and NFE GnomAD figures (including and excluding variants with AF > 0.005), as well as allele freq ratios

ratios <- function(vcf){
  mutate(vcf, 
         "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,
         "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_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_FHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats)

masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_noFHx_goodQ_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats)

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats)

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios <- 
  ratios(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats)

Create data frames with Agilent SureSelect whole exome genes (all ENST, protein-coding only and non-protein-coding) for volcano 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_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_list <- tibble("SYMBOL"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$SYMBOL,"Gene"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$Gene,"BIOTYPE"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)
allGenes_masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_list <- tibble("SYMBOL"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$SYMBOL,"Gene"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$Gene,"BIOTYPE"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  arrange(SYMBOL)
allGenes_masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype2_list <- tibble("SYMBOL"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$SYMBOL,"Gene"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$Gene,"BIOTYPE"=masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  filter(!BIOTYPE%in%c("protein_coding")) %>% 
  arrange(SYMBOL)

volcano <- function(vcf,gene_list,GnomAD,Total_Sample_Alleles){
  water <- vcf %>% 
  mutate(Total_Case_Alleles=Total_Sample_Alleles) %>% 
  select(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) %>% 
  distinct(SYMBOL,Gene,Total_Gene_Count,Total_Case_Alleles) 
  water2 <- bind_rows(water, anti_join(gene_list,water,by="Gene")) %>% 
  mutate_at(vars("Total_Gene_Count"),funs(replace(., is.na(.), 0))) %>% 
  mutate_at(vars("Total_Case_Alleles"),funs(replace(., is.na(.), Total_Sample_Alleles)))
  water3 <- right_join(water2,GnomAD,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)
  return(water3)
}

volcano_FHx_sampleAF0.01_allGenes <- volcano(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all,Total_Sample_Alleles_FHx)

volcano_FHx_sampleAF0.01_biotype_allGenes <- volcano(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list_protein_coding_only,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding,Total_Sample_Alleles_FHx)

volcano_FHx_sampleAF0.01_biotype2_allGenes <- volcano(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list_non_protein_coding_only,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding,Total_Sample_Alleles_FHx)

volcano_noFHx_sampleAF0.01_allGenes <- volcano(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_all,Total_Sample_Alleles_noFHx)

volcano_noFHx_sampleAF0.01_biotype_allGenes <- volcano(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list_protein_coding_only,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_protein_coding,Total_Sample_Alleles_noFHx)

volcano_noFHx_sampleAF0.01_biotype2_allGenes <- volcano(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,allGenes_AgilentSSv6_list_non_protein_coding_only,GnomAD_stats_AgilentSSv6_LOF_VEP_ENSTcanonical_non_protein_coding,Total_Sample_Alleles_noFHx)

Calculate Fisher’s test values (Ensembl Canonical)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,Total_Case_Alleles=Total_Sample_Alleles_FHx) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_FHx_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_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_ENST <- apply(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,230,244,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_FHx_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_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005 = sapply(fisherresults_ENST,function(x) round(x$estimate,4))
masterfile_eoc_FHx_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="~"))
masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$BH_0.005 = p.adjust(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005,method = "BH")

fisherresults2_ENST <- apply(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,233,244,246)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_FHx_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_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005_NFE = sapply(fisherresults2_ENST,function(x) round(x$estimate,4))
masterfile_eoc_FHx_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="~"))
masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$BH_0.005_NFE = p.adjust(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005_NFE,method = "BH")

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios,Total_Case_Alleles=Total_Sample_Alleles_noFHx) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_noFHx_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_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_ENST <- apply(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,230,244,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_noFHx_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_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005 = sapply(fisherresults_ENST,function(x) round(x$estimate,4))
masterfile_eoc_noFHx_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="~"))
masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$BH_0.005 = p.adjust(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005,method = "BH")

fisherresults2_ENST <- apply(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,233,244,246)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

masterfile_eoc_noFHx_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_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$OR_0.005_NFE = sapply(fisherresults2_ENST,function(x) round(x$estimate,4))
masterfile_eoc_noFHx_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="~"))
masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$BH_0.005_NFE = p.adjust(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2$P_value_Fisher_0.005_NFE,method = "BH")

Calculate Fisher’s test values (RefSeq Canonical)

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios,Total_Case_Alleles=Total_Sample_Alleles_FHx) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_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_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_NM <- apply(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,230,244,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

fisherresults2_NM <- apply(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,233,244,246)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios,Total_Case_Alleles=Total_Sample_Alleles_noFHx) %>% 
  left_join(GnomAD_stats_LOF_VEP[,c(2,6,10,14:29)],by="Feature",copy=FALSE) 
masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2 <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults, MAX_AN=replace(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN, is.na(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN), max(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_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_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults$MAX_AN_NFE, na.rm=TRUE)))

fisherresults_NM <- apply(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,230,244,245)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

fisherresults2_NM <- apply(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults2[,c(222,233,244,246)],1,function(x) {fisher.test(rbind(x[1:2],c(x[3]-x[1],x[4]-x[2])))})

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

Calculate Fisher’s test values (volcano plot)

fisherresults_allGenes_volcano <- apply(volcano_FHx_sampleAF0.01_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_allGenes_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_allGenes$OR_0.005 = sapply(fisherresults_allGenes_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_allGenes$`95%CI_0.005` = sapply(fisherresults_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_allGenes$BH_0.005 = p.adjust(volcano_FHx_sampleAF0.01_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_NFE_volcano <- apply(volcano_FHx_sampleAF0.01_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_NFE_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_allGenes$OR_0.005_NFE = sapply(fisherresults_NFE_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_allGenes$BH_0.005_NFE = p.adjust(volcano_FHx_sampleAF0.01_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

fisherresults_biotype_allGenes_volcano <- apply(volcano_FHx_sampleAF0.01_biotype_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_biotype_allGenes_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_biotype_allGenes$OR_0.005 = sapply(fisherresults_biotype_allGenes_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_biotype_allGenes$`95%CI_0.005` = sapply(fisherresults_biotype_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_biotype_allGenes$BH_0.005 = p.adjust(volcano_FHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_biotype_NFE_volcano <- apply(volcano_FHx_sampleAF0.01_biotype_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_biotype_NFE_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_biotype_allGenes$OR_0.005_NFE = sapply(fisherresults_biotype_NFE_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_biotype_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_biotype_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_biotype_allGenes$BH_0.005_NFE = p.adjust(volcano_FHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

fisherresults_biotype2_allGenes_volcano <- apply(volcano_FHx_sampleAF0.01_biotype2_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_biotype2_allGenes_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_biotype2_allGenes$OR_0.005 = sapply(fisherresults_biotype2_allGenes_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_biotype2_allGenes$`95%CI_0.005` = sapply(fisherresults_biotype2_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_biotype2_allGenes$BH_0.005 = p.adjust(volcano_FHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_biotype2_NFE_volcano <- apply(volcano_FHx_sampleAF0.01_biotype2_allGenes[,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]))))})

volcano_FHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_biotype2_NFE_volcano,function(x) round(x$p.value,15))
volcano_FHx_sampleAF0.01_biotype2_allGenes$OR_0.005_NFE = sapply(fisherresults_biotype2_NFE_volcano,function(x) round(x$estimate,15))
volcano_FHx_sampleAF0.01_biotype2_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_biotype2_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_FHx_sampleAF0.01_biotype2_allGenes$BH_0.005_NFE = p.adjust(volcano_FHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

fisherresults_allGenes_volcano <- apply(volcano_noFHx_sampleAF0.01_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_allGenes_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_allGenes$OR_0.005 = sapply(fisherresults_allGenes_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_allGenes$`95%CI_0.005` = sapply(fisherresults_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_allGenes$BH_0.005 = p.adjust(volcano_noFHx_sampleAF0.01_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_NFE_volcano <- apply(volcano_noFHx_sampleAF0.01_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_NFE_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_allGenes$OR_0.005_NFE = sapply(fisherresults_NFE_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_allGenes$BH_0.005_NFE = p.adjust(volcano_noFHx_sampleAF0.01_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

fisherresults_biotype_allGenes_volcano <- apply(volcano_noFHx_sampleAF0.01_biotype_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_biotype_allGenes_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_biotype_allGenes$OR_0.005 = sapply(fisherresults_biotype_allGenes_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_biotype_allGenes$`95%CI_0.005` = sapply(fisherresults_biotype_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_biotype_allGenes$BH_0.005 = p.adjust(volcano_noFHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_biotype_NFE_volcano <- apply(volcano_noFHx_sampleAF0.01_biotype_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_biotype_NFE_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_biotype_allGenes$OR_0.005_NFE = sapply(fisherresults_biotype_NFE_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_biotype_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_biotype_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_biotype_allGenes$BH_0.005_NFE = p.adjust(volcano_noFHx_sampleAF0.01_biotype_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

fisherresults_biotype2_allGenes_volcano <- apply(volcano_noFHx_sampleAF0.01_biotype2_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005 = sapply(fisherresults_biotype2_allGenes_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$OR_0.005 = sapply(fisherresults_biotype2_allGenes_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$`95%CI_0.005` = sapply(fisherresults_biotype2_allGenes_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$BH_0.005 = p.adjust(volcano_noFHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005,method = "BH")

fisherresults_biotype2_NFE_volcano <- apply(volcano_noFHx_sampleAF0.01_biotype2_allGenes[,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]))))})

volcano_noFHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005_NFE = sapply(fisherresults_biotype2_NFE_volcano,function(x) round(x$p.value,15))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$OR_0.005_NFE = sapply(fisherresults_biotype2_NFE_volcano,function(x) round(x$estimate,15))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$`95%CI_0.005_NFE` = sapply(fisherresults_biotype2_NFE_volcano,function(x) paste(round(x$conf.int[1:2],2),collapse="~"))
volcano_noFHx_sampleAF0.01_biotype2_allGenes$BH_0.005_NFE = p.adjust(volcano_noFHx_sampleAF0.01_biotype2_allGenes$P_value_Fisher_0.005_NFE,method = "BH")

Identify variants with filtering AF > maximum credible population AF in GnomAD for hereditary EOC (0.000181)

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_FHx_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))

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_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_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_noFHx_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))

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults3 <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_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))

Identify variants that PASS/FAIL the LOFTEE 50bp rule

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4 <- mutate(masterfile_eoc_FHx_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)))

masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4 <- mutate(masterfile_eoc_FHx_goodQ2_HIGH_NMcanonical_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)))

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4 <- mutate(masterfile_eoc_noFHx_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)))

masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_sampleAF0.01_withGnomADstats_ratios_fisherresults4 <- mutate(masterfile_eoc_noFHx_goodQ2_HIGH_NMcanonical_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 and annotate variants accordingly

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

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

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

Remove redundant objects and output data

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

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

select(masterfile_eoc_FHx_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="FHx_HIGH_ENSTcanonical_BH_p_values.tsv")

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

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

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

select(masterfile_eoc_noFHx_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="FHx_HIGH_ENSTcanonical_BH_p_values.tsv")

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

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

Calculate total LoF variants in sample vs total LoF variants in GnomAD non-cancer NFE and use for chi-squared test

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("LoF", "No LoF", "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)
}

a <- sum(volcano_FHx_sampleAF0.01_allGenes$Total_Gene_Count) %>% as.numeric()
b <- Total_Sample_Alleles_FHx*(nrow(allGenes_AgilentSSv6_list)) %>% as.numeric()
b1 <- b-a %>% as.numeric()
c <- sum(volcano_FHx_sampleAF0.01_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
d <- (max(volcano_FHx_sampleAF0.01_allGenes$MAX_AN))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
d1 <- d-c %>% as.numeric()
e <- sum(volcano_FHx_sampleAF0.01_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
f <- (max(volcano_FHx_sampleAF0.01_allGenes$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
f1 <- f-e %>% as.numeric()

g <- sum(volcano_FHx_sampleAF0.01_biotype_allGenes$Total_Gene_Count) %>% as.numeric()
h <- Total_Sample_Alleles_FHx*(nrow(allGenes_AgilentSSv6_list_protein_coding_only)) %>% as.numeric()
h1 <- h-g %>% as.numeric()
i <- sum(volcano_FHx_sampleAF0.01_biotype_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
j <- (max(volcano_FHx_sampleAF0.01_biotype_allGenes$MAX_AN))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
j1 <- j-i %>% as.numeric()
k <- sum(volcano_FHx_sampleAF0.01_biotype_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
l <- (max(volcano_FHx_sampleAF0.01_biotype_allGenes$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
l1 <- l-k %>% as.numeric()

chiX2_0.005_FHx <- matrix(c(a,b1,c,d1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_FHx <- matrix(c(a,b1,e,f1),nrow=2,byrow=TRUE)
chiX2_0.005_biotype_FHx <- matrix(c(g,h1,i,j1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_biotype_FHx <- matrix(c(g,h1,k,l1),nrow=2,byrow=TRUE)

m <- sum(volcano_noFHx_sampleAF0.01_allGenes$Total_Gene_Count) %>% as.numeric()
n <- Total_Sample_Alleles_noFHx*(nrow(allGenes_AgilentSSv6_list)) %>% as.numeric()
n1 <- n-m %>% as.numeric()
o <- sum(volcano_noFHx_sampleAF0.01_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
p <- (max(volcano_noFHx_sampleAF0.01_allGenes$MAX_AN))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
p1 <- p-o %>% as.numeric()
q <- sum(volcano_noFHx_sampleAF0.01_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
r <- (max(volcano_noFHx_sampleAF0.01_allGenes$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list) %>% as.numeric()
r1 <- r-q %>% as.numeric()

s <- sum(volcano_noFHx_sampleAF0.01_biotype_allGenes$Total_Gene_Count) %>% as.numeric()
t <- Total_Sample_Alleles_noFHx*(nrow(allGenes_AgilentSSv6_list_protein_coding_only)) %>% as.numeric()
t1 <- t-s %>% as.numeric()
u <- sum(volcano_noFHx_sampleAF0.01_biotype_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005) %>% as.numeric()
v <- (max(volcano_noFHx_sampleAF0.01_biotype_allGenes$MAX_AN))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
v1 <- v-u %>% as.numeric()
w <- sum(volcano_noFHx_sampleAF0.01_biotype_allGenes$FILTER_RF_LOF_HIGHIMPACT_AC_0.005_NFE) %>% as.numeric()
x <- (max(volcano_noFHx_sampleAF0.01_biotype_allGenes$MAX_AN_NFE))*nrow(allGenes_AgilentSSv6_list_protein_coding_only) %>% as.numeric()
x1 <- x-w %>% as.numeric()

chiX2_0.005_noFHx <- matrix(c(m,n1,o,p1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_noFHx <- matrix(c(m,n1,q,r1),nrow=2,byrow=TRUE)
chiX2_0.005_biotype_noFHx <- matrix(c(s,t1,u,v1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_biotype_noFHx <- matrix(c(s,t1,w,x1),nrow=2,byrow=TRUE)

chiX2_0.005_NFE_FHx_noFHx <- matrix(c(chiX2_0.005_NFE_FHx[1,],chiX2_0.005_NFE_noFHx[1,]),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_biotype_FHx_noFHx <- matrix(c(chiX2_0.005_NFE_biotype_FHx[1,],chiX2_0.005_NFE_biotype_noFHx[1,]),nrow=2,byrow=TRUE)

setwd("~/OneDrive - The University of Melbourne/R data")
sink(file="masterfile_eoc_LOF_volcano_chiX2_results.txt",append=TRUE)

chisq.test(chiX2_0.005_FHx,correct=FALSE)
oddsratio(chiX2_0.005_FHx)
chisq.test(chiX2_0.005_NFE_FHx,correct=FALSE)
oddsratio(chiX2_0.005_NFE_FHx)
chisq.test(chiX2_0.005_biotype_FHx,correct=FALSE)
oddsratio(chiX2_0.005_biotype_FHx)
chisq.test(chiX2_0.005_NFE_biotype_FHx,correct=FALSE)
oddsratio(chiX2_0.005_NFE_biotype_FHx)

chisq.test(chiX2_0.005_noFHx,correct=FALSE)
oddsratio(chiX2_0.005_noFHx)
chisq.test(chiX2_0.005_NFE_noFHx,correct=FALSE)
oddsratio(chiX2_0.005_NFE_noFHx)
chisq.test(chiX2_0.005_biotype_noFHx,correct=FALSE)
oddsratio(chiX2_0.005_biotype_noFHx)
chisq.test(chiX2_0.005_NFE_biotype_noFHx,correct=FALSE)
oddsratio(chiX2_0.005_NFE_biotype_noFHx)

chisq.test(chiX2_0.005_NFE_FHx_noFHx,correct=FALSE)
chisq.test(chiX2_0.005_NFE_biotype_FHx_noFHx,correct=FALSE)

sink()

Remove variants with non-protein-coding biotypes

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

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

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

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

Calculate number of retained protein-coding LoF variants per individual in FHx/noFHx cohorts, and output results.

masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- select(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype,Sample,CHROM,POS,REF,ALT,Sample.GT,SYMBOL,Gene,Feature,HGVSc)
samples_FHx <- masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson[,1] %>% unique()
samples_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- data.frame(Sample=factor(),No_of_Variants=double())
for(sample in samples_FHx){
  sample_data <- filter(masterfile_eoc_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,Sample%in%c(sample))
  sample_variants <- variant_counts_AFs(sample_data,Total_Sample_Alleles_FHx)
  no_of_variants <- sum(sample_variants$Total_Allele_Count)
  samples_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- add_row(samples_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,Sample=sample,No_of_Variants=no_of_variants)
}
write_tsv(samples_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,path="samples_FHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson.tsv")

masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- select(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype,Sample,CHROM,POS,REF,ALT,Sample.GT,SYMBOL,Gene,Feature,HGVSc)
samples_noFHx <- masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson[,1] %>% unique()
samples_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- data.frame(Sample=factor(),No_of_Variants=double())
for(sample in samples_noFHx){
  sample_data <- filter(masterfile_eoc_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,Sample%in%c(sample))
  sample_variants <- variant_counts_AFs(sample_data,Total_Sample_Alleles_noFHx)
  no_of_variants <- sum(sample_variants$Total_Allele_Count)
  samples_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson <- add_row(samples_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,Sample=sample,No_of_Variants=no_of_variants)
}
write_tsv(samples_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson,path="samples_noFHx_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_variantsPerPerson.tsv")

Output data for log p-value volcano plots

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

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

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

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

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

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