WES VAF Analysis: Patient1 and Cell Lines L1-L2

1 Load libraries

library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
library(viridis)

# Optional: theme for plots
theme_set(theme_bw())

2 Load AD files

# Define files

files <- list(
N1 = "N1.AD.txt",
T1 = "T1.AD.txt",
X1 = "X1.AD.txt",
XX1 = "XX1.AD.txt"
)

# Function to read and process AD file

read_ad <- function(file){
read_tsv(file, col_names = c("CHROM","POS","REF","ALT","AD")) %>%
separate(AD, into=c("ref_count","alt_count"), sep=",", convert=TRUE) %>%
mutate(VAF = alt_count / (ref_count + alt_count))
}

# Load all samples into a list

ad_list <- lapply(files, read_ad)

# Name list elements

names(ad_list) <- names(files)

3 Compute basic statistics

# Compute basic statistics

vaf_summary <- lapply(ad_list, function(df){
summary(df$VAF)
})

vaf_summary

$N1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.3596  0.4902  0.5537  0.9951  1.0000       1 

$T1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.3529  0.5000  0.5663  0.9942  1.0000 

$X1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.3626  0.5000  0.5692  1.0000  1.0000       1 

$XX1
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.0000  0.3448  0.5000  0.5650  1.0000  1.0000       3 

4 Plot VAF distributions

# Combine all VAFs for plotting

vaf_df <- bind_rows(
lapply(names(ad_list), function(n) {
ad_list[[n]] %>% select(VAF) %>% mutate(Sample = n)
})
)

ggplot(vaf_df, aes(x=VAF, fill=Sample)) +
geom_histogram(bins=50, alpha=0.6, position="identity") +
scale_fill_viridis(discrete = TRUE) +
labs(title="VAF Distribution Across Samples",
x="Variant Allele Frequency (VAF)",
y="Number of Variants") +
theme_minimal()

NA
NA
NA

5 Identify shared variants

library(dplyr)

# Remove duplicate variants by CHROM and POS
T1_unique <- ad_list$T1 %>% distinct(CHROM, POS, .keep_all = TRUE)
XX1_unique <- ad_list$XX1 %>% distinct(CHROM, POS, .keep_all = TRUE)

# Find shared variants between tumor (T1) and cell line (XX1)
shared_T1_XX1 <- inner_join(
  T1_unique %>% select(CHROM, POS),
  XX1_unique %>% select(CHROM, POS),
  by = c("CHROM", "POS"),
  relationship = "many-to-many"  # explicitly allow many-to-many
)

cat("Number of shared variants between T1 and XX1:", nrow(shared_T1_XX1), "\n")

Number of shared variants between T1 and XX1: 44650 

6 Compute shared vs unique variants

get_variant_id <- function(df){
  df %>% mutate(VariantID = paste(CHROM, POS, sep="_")) %>% select(VariantID)
}

variant_ids <- lapply(ad_list, get_variant_id)
variant_ids <- lapply(variant_ids, pull)

shared_T1_X1 <- intersect(variant_ids$T1, variant_ids$X1)
shared_T1_XX1 <- intersect(variant_ids$T1, variant_ids$XX1)

unique_T1 <- setdiff(variant_ids$T1, union(variant_ids$X1, variant_ids$XX1))
unique_X1 <- setdiff(variant_ids$X1, variant_ids$T1)
unique_XX1 <- setdiff(variant_ids$XX1, variant_ids$T1)

variant_summary <- data.frame(
  Comparison = c("Shared_T1_X1","Shared_T1_XX1","Unique_T1","Unique_X1","Unique_XX1"),
  Count = c(length(shared_T1_X1), length(shared_T1_XX1),
            length(unique_T1), length(unique_X1), length(unique_XX1))
)
variant_summary

library(VennDiagram)
library(grid)

variants_sets <- list(
  T1 = variant_ids$T1,
  X1 = variant_ids$X1,
  XX1 = variant_ids$XX1
)

venn.diagram(
  x = variants_sets,
  filename = "Venn_Tumor_CellLines.png",
  fill = c("red", "green", "blue"),
  alpha = 0.5,
  cex = 1.5,
  cat.cex = 1.2,
  main = "Overlap of Variants: Tumor and Cell Lines"
)

[1] 1

# Use grid.draw to render
grid.newpage()
grid.draw(venn.plot)

7 VAF Heatmap for Shared Variants

library(tidyverse)
library(pheatmap)

# Combine all AD data
ad_df <- bind_rows(
  lapply(names(ad_list), function(s) {
    ad_list[[s]] %>% mutate(Sample = s)
  })
) %>%
  mutate(VariantID = paste(CHROM, POS, sep="_"))

# Shared variants between tumor and at least one cell line
shared_variants <- union(shared_T1_X1, shared_T1_XX1)

# FIX: Remove duplicates before pivoting
vaf_df <- ad_df %>%
  filter(VariantID %in% shared_variants) %>%
  group_by(VariantID, Sample) %>%      
  summarise(VAF = mean(VAF, na.rm = TRUE), .groups = "drop") %>%  # FIX
  pivot_wider(names_from = Sample, values_from = VAF, values_fill = 0)

# Convert to matrix
vaf_matrix <- as.matrix(vaf_df[,-1])
rownames(vaf_matrix) <- vaf_df$VariantID

# Heatmap
pheatmap(
  vaf_matrix,
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  color = colorRampPalette(c("white","red"))(50),
  main = "VAF Heatmap: Shared Variants Tumor vs Cell Lines"
)

NA
NA
NA