Thesis Synonymous Variants 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 with GnomAD non-cancer AF > 0.005

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

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

sample_ethnicities <- read.delim("exomes_pca_with_gnomad_info.tsv", stringsAsFactors=FALSE)
NFE_samples <- filter(sample_ethnicities, GnomAD_class%in%"European")

ViP_Complete_Cohort <- read.delim("ViP_LoF_Complete_Cohort.txt", stringsAsFactors=FALSE)
ViP_Complete_Cohort_list_NFE <- filter(ViP_Complete_Cohort,Exome.ID%in%c(NFE_samples$Sample))

masterfile_eoc <- filter(masterfile_rare_non_cancer,Sample%in%c(ViP_Complete_Cohort_list_NFE$Exome.ID)) %>% 
  filter(Sample%in%(NFE_samples$Sample)) %>%
  filter(GnomAD_v2.1_non_cancer_AF_nfe<=0.005)

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)
rm(masterfile_eoc)

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

masterfile_eoc_goodQ <- filter(masterfile_eoc_withPath, (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))

rm(masterfile_eoc_withPath)

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_goodQ2 <- filter(masterfile_eoc_goodQ,Sample%in%c(ViP_Discovery_Cohort_list))

Exclude HIGH and MODERATE impact variants

masterfile_eoc_goodQ2_LOW <- filter(masterfile_eoc_goodQ2, IMPACT%in%"LOW")

rm(masterfile_eoc_goodQ2)

Exclude non-synonymous variants

masterfile_eoc_goodQ2_LOW_synonymous_NoLOFTEE <- filter(masterfile_eoc_goodQ2_LOW,Consequence%in%c("synonymous_variant","synonymous_variant,NMD_transcript_variant")) %>% 
  filter(!LoF%in%c("HC","LC"))

rm(masterfile_eoc_goodQ2_LOW)

Separate Ensembl canonical and RefSeq canonical variants

masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical <- filter(masterfile_eoc_goodQ2_LOW_synonymous_NoLOFTEE,CANONICAL%in%"YES") 

Sample variant counts and AFs

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

Exclude variants with sample AF >0.01

masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01 <- filter(masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_2,Sample_AF < 0.01)

Output list of genes and variants with sample AF >0.01

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

Sample gene counts and frequencies

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

Add GnomAD gene-level data

ensembl_biotypes <- read.delim("ensembl_biotypes.tsv", stringsAsFactors=FALSE) %>% 
  select(ensembl_transcript_id,transcript_biotype) %>% 
  dplyr::rename("Feature"="ensembl_transcript_id","BIOTYPE"="transcript_biotype")

all_features_exons_only_fraction_covered_agilent_v6_base_figures <- read.delim("~/all_features_exons_only_fraction_covered_agilent_v6_150bp_extended.tsv")

gnomad_coverage_coding_features <- read.delim("~/gnomad_coverage_coding_features.tsv") %>% 
  mutate("total_10x_non_cancer_count"=((exome_10x)*56885)+((genome_10x)*7718)) %>% 
  mutate("total_10x_non_cancer_NFE_count"=((exome_10x)*51377)+((genome_10x)*7718))
gnomad_coverage_coding_features$total_10x_non_cancer_count=round(gnomad_coverage_coding_features$total_10x_non_cancer_count)
gnomad_coverage_coding_features$total_10x_non_cancer_NFE_count=round(gnomad_coverage_coding_features$total_10x_non_cancer_NFE_count)

GnomAD_stats_synonymous_VEP_NO_RF <- read.delim("agilent_v6_gnomad_v2.1.1_genestats_canonical_SYN_noLOFTEE.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  left_join(ensembl_biotypes,by="Feature",copy=FALSE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  distinct(.keep_all = FALSE) %>% 
  left_join(all_features_exons_only_fraction_covered_agilent_v6_base_figures[,c(1,4)],by="Feature",copy=FALSE) %>% 
  left_join(gnomad_coverage_coding_features[,c(1,6,7)],by="Feature",copy=FALSE)
  
masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats <- left_join(masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_2, GnomAD_stats_synonymous_VEP_NO_RF[,c(2,20:25,27:29)],by="Feature",copy=FALSE)
masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats2 <- mutate_at(masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats,vars(starts_with("FILTER_")),funs(replace(., is.na(.), 0))) %>% 
  mutate_if(grepl("popmax$", names(.)),funs(ifelse(. == "NA", 0, as.numeric(.)))) %>% 
  mutate_at(vars(ends_with("popmax")),funs(replace(., is.na(.), 0)))

Create data frames with Agilent SureSelect whole exome genes (all ENST, protein-coding only and non-protein-coding)

GnomAD_stats_AgilentSSv6_list <-  read.delim("agilent_v6_gnomad_v2.1.1_genestats_canonical_SYN_noLOFTEE.tsv", header=TRUE, row.names=NULL, stringsAsFactors=FALSE) %>% 
  filter(CANONICAL=="YES") %>% 
  distinct(.keep_all = FALSE) %>% 
  select(SYMBOL:ENSP)

GnomAD_stats_AgilentSSv6_SYN_NO_LOFTEE_VEP_ENSTcanonical_protein_coding <- filter(GnomAD_stats_synonymous_VEP_NO_RF,Feature%in%c(GnomAD_stats_AgilentSSv6_list$Feature)) %>% 
  dplyr::rename("Gene"="ENSG") %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  distinct(.keep_all = FALSE)

allGenes_AgilentSSv6_list_protein_coding_only <- tibble("SYMBOL"=GnomAD_stats_AgilentSSv6_SYN_NO_LOFTEE_VEP_ENSTcanonical_protein_coding$SYMBOL,"Gene"=GnomAD_stats_AgilentSSv6_SYN_NO_LOFTEE_VEP_ENSTcanonical_protein_coding$Gene,"Feature"=GnomAD_stats_AgilentSSv6_SYN_NO_LOFTEE_VEP_ENSTcanonical_protein_coding$Feature) %>% 
  distinct(.keep_all=TRUE) %>% 
  arrange(SYMBOL)

allGenes_masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_biotype_list <- tibble("SYMBOL"=masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats2$SYMBOL,"Gene"=masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats2$Gene,"BIOTYPE"=masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats2$BIOTYPE) %>% 
  distinct(.keep_all=TRUE) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  arrange(SYMBOL)

comparison_sampleAF0.01_biotype_allGenes <- masterfile_eoc_goodQ_0.005_LOW_synonymous_NoLOFTEE_ENSTcanonical_sampleAF0.01_withGnomADstats2 %>% 
  mutate(Total_Case_Bases=(n_distinct(masterfile_eoc_goodQ2_LOW_synonymous_NoLOFTEE$Sample)*2*Total_Bases_Covered)) %>% 
  filter(BIOTYPE%in%c("protein_coding")) %>% 
  select(SYMBOL,Gene,Total_Gene_Count,Total_Case_Bases) %>% 
  distinct(SYMBOL,Gene,Total_Gene_Count,Total_Case_Bases) 
comparison_sampleAF0.01_biotype_allGenes2 <- bind_rows(comparison_sampleAF0.01_biotype_allGenes, anti_join(allGenes_AgilentSSv6_list_protein_coding_only,comparison_sampleAF0.01_biotype_allGenes,by="Gene")) %>% 
  mutate_at(vars("Total_Gene_Count"),funs(replace(., is.na(.), 0)))
comparison_sampleAF0.01_biotype_allGenes3 <- right_join(
  comparison_sampleAF0.01_biotype_allGenes2,
  GnomAD_stats_AgilentSSv6_SYN_NO_LOFTEE_VEP_ENSTcanonical_protein_coding,
  by="Gene",copy=FALSE) %>% 
  mutate_at(vars("Total_Case_Bases"),funs(replace(., is.na(.), n_distinct(masterfile_eoc_goodQ2_LOW_synonymous_NoLOFTEE$Sample)*2*Total_Bases_Covered))) %>% 
  select("SYMBOL.x","Gene","Total_Gene_Count","Total_Case_Bases","FILTER_RF_SYN_NO_LOFTEE_AC_0.005","FILTER_RF_SYN_NO_LOFTEE_AC_0.005_NFE","Total_Bases_Covered","total_10x_non_cancer_count","total_10x_non_cancer_NFE_count") %>% 
  dplyr::rename("SYMBOL"="SYMBOL.x") %>% 
  mutate("Total_GnomAD_Bases"=total_10x_non_cancer_count*2*Total_Bases_Covered) %>% 
  mutate("Total_GnomAD_NFE_Bases"=total_10x_non_cancer_NFE_count*2*Total_Bases_Covered) %>% 
  filter(Total_Bases_Covered!=0) %>% 
  distinct(.keep_all=TRUE)

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("Syn", "No Syn", "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(comparison_sampleAF0.01_biotype_allGenes3$Total_Gene_Count) %>% as.numeric()
b <- n_distinct(masterfile_eoc_goodQ2_LOW_synonymous_NoLOFTEE$Sample)*2*33218304 %>% as.numeric()
b1 <- b-a %>% as.numeric()
c <- sum(comparison_sampleAF0.01_biotype_allGenes3$FILTER_RF_SYN_NO_LOFTEE_AC_0.005) %>% as.numeric()
d <- (30651075*2*118479)+(32482680*2*15708) %>% as.numeric()
d1 <- d-c %>% as.numeric()
e <- sum(comparison_sampleAF0.01_biotype_allGenes3$FILTER_RF_SYN_NO_LOFTEE_AC_0.005_NFE) %>% as.numeric()
f <- (30651075*2*51377)+(32482680*2*7718) %>% as.numeric()
f1 <- f-e %>% as.numeric()

chiX2_0.005_biotype <- matrix(c(a,b1,c,d1),nrow=2,byrow=TRUE)
chiX2_0.005_NFE_biotype <- matrix(c(a,b1,e,f1),nrow=2,byrow=TRUE)

sink(file="masterfile_eoc_SYN_NO_LOFTEE_comparison_chiX2_results.txt",append=FALSE)

chisq.test(chiX2_0.005_biotype,correct=FALSE)
oddsratio(chiX2_0.005_biotype)
chisq.test(chiX2_0.005_NFE_biotype,correct=FALSE)
oddsratio(chiX2_0.005_NFE_biotype)

sink()