Thesis LoF Variant Burden Analysis Script

THIS SCRIPT REQUIRES THE LOF VARIANT FILTERING AND RANKING SCRIPT (https://rpubs.com/deepsubs/nature_comms_paper_2020) TO BE RUN FIRST, AND THE SPECIFIED OUTPUT FILE FROM THAT SCRIPT TO BE SAVED IN THE SAME DIRECTORY AS THIS SCRIPT

Load tidyverse package

library("tidyverse")

Load data from LoF variant filtering and ranking analysis

masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01 <- read.csv("masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01.csv")

Remove variants with non-protein-coding biotypes

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

Remove variants with GnomAD PopMax AF > 0.005

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

Remove variants with ratio < 3

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

Output values for binomial test

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_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson <- select(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3,Sample,CHROM,POS,REF,ALT,Sample.GT,SYMBOL,Gene,Feature,HGVSc)
samples <- masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson[,1] %>% unique()
samples_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson <- data.frame(Sample=factor(),No_of_Variants=double())
for(sample in samples){
  sample_data <- filter(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson,Sample%in%c(sample))
  sample_variants <- variant_counts_AFs(sample_data,(n_distinct(samples)*2))
  no_of_variants <- sum(sample_variants$Total_Allele_Count)
  samples_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson <- add_row(samples_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson,Sample=sample,No_of_Variants=no_of_variants)
}

binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3 <- samples_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3_variantsPerPerson %>%
  rename(gene_Var = "No_of_Variants") %>%
  count(gene_Var) %>%
  rename(obs_n_samples = "n") %>%
  add_row(gene_Var = 0, obs_n_samples = (510)-n_distinct(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3$Sample)) %>%
  arrange(gene_Var) %>%
  mutate(Total_No=gene_Var*obs_n_samples) %>%
  mutate(p_LoF_variant_perSample=sum(Total_No)/(510*2*(n_distinct(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3$Gene)))) %>%
  mutate(binomial_distrib_p=(dbinom(gene_Var,2*(n_distinct(masterfile_eoc_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3$Gene)),p_LoF_variant_perSample,log=FALSE))) %>%
  mutate(exp_n_samples = round(binomial_distrib_p*(510),digits=0)) %>%
  write_tsv(path="binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3.tsv")

multinomial <- binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3[c(1:8),c(2,5,6)] %>%
  add_row(obs_n_samples=(sum(binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3[c(9:13),2])),
          binomial_distrib_p=(sum(binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3[c(9:13),5])),
          exp_n_samples=(sum(binomial_results_goodQ2_HIGH_ENSTcanonical_sampleAF0.01_biotype_GnomADpopmax0.005_ratio3[c(9:13),6])))

Perform chi-squared test

chisq.test(x=multinomial$obs_n_samples,
           correct=FALSE,
           p=multinomial$binomial_distrib_p,
           rescale.p=TRUE)

    Chi-squared test for given probabilities

data:  multinomial$obs_n_samples
X-squared = 10.56, df = 8, p-value = 0.2279