CNV_analysis.Rmd

---
title: "Extended Data Figure 3b CNV analysis"
author: "Sam Buckberry"
date: "`r Sys.Date()`"
output:
    html_document:
        code_folding: show
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, message=FALSE}
source("R/project_functions.R")
```


## Experimental aim

The aim of this experiment is to determine if naive treatment of cells during reprogramming causes any genomic instability, compared with conventional primed reprogramming. 

We performed Primed, TNT and NtP reprogramming for two independent fibroblast lines (32F and 38F). To test for genomic deletions or insertions, we extracted DNA and performed genome-wide genotyping using the Illumina Global Screening array v2 with 663,320 SNPs surveyed. Data were processed using GSGT Version 1.9.4. 

```{r, echo=FALSE, cache=TRUE}
# Add sample IDs
ids <- c(TNT_32F_P3="19C106621_KY1_MA1764-68_FinalReport.txt.gz",
  TNT_38F_P3="19C106622_KY2_MA1764-72_FinalReport.txt.gz",
  NtP_38F_P12_P16="19C106623_KY3_MA1764-78_FinalReport.txt.gz",
  Primed_32F_P33="19C108849_32FP33_MA1796-78_FinalReport.txt.gz",
  Primed_38F_P22="19C108850_38FP22_MA1796-79_FinalReport.txt.gz",
  NtP_32F_P17_P23="19C108851_32FP1723_MA1796-85_FinalReport.txt.gz")

sample_info <- data.frame(ID=names(ids), File=ids)

data.table(sample_info)
```

## Data analysis

List the data files and check ids. 
```{r, cache=TRUE}
snp_fls <- list.files(path = "snp_array/",
                      pattern = "FinalReport.txt.gz",
                      full.names = TRUE)


if (all(basename(snp_fls) == ids)){
    names(snp_fls) <- names(ids)
} else {
    stop()
}
```

Function to read snp chip data and format for analysis
```{r, cache=TRUE}
read_snp <- function(path, chroms=c(1:22, "X")){
    
    df <- fread(input = path, skip = 10, header = TRUE)
    
    colnames(df) <- str_replace_all(colnames(df),
                                    pattern = " ",
                                    replacement = "_")
    df <- df[df$Chr %in% chroms, ]
    
    ind <- match(basename(path), ids)
    
    df$Sample_ID <- names(ids)[ind]
    
    return(df)
}
```


Function to plot allele frequencies for each chromosome
```{r, cache=TRUE}

plot_BAF <- function(snp_fls, chrom_list=chroms){
    
    df <- lapply(snp_fls, read_snp, chroms=chrom_list) %>% do.call(rbind, .)
    df$Chr <- factor(df$Chr, levels=c(1:22, "X", "Y"))
    
    gg_freq <- ggplot(data = df, aes(y = B_Allele_Freq, x = Position)) +
        geom_point(size=0.1, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0)) +
        xlab("Position") + ylab("BAF") +
        facet_grid(Sample_ID~Chr, scales = "free", space = "free") +
        sams_pub_theme() +
        theme(panel.spacing.x = unit(0, "lines")) +
        theme(panel.border = element_rect(fill = NA))
    
    return(gg_freq)
}

```

Function to plot Log Ratio of alleles
```{r, cache=TRUE}
plot_LR <- function(snp_fls, chrom_list=chroms){
    
    df <- lapply(snp_fls, read_snp, chroms=chrom_list) %>% do.call(rbind, .)
    df$Chr <- factor(df$Chr, levels=c(1:22, "X", "Y"))
    
    gg_log_ratio <- ggplot(df, aes(x = Position, y = Log_R_Ratio)) +
        geom_point(size=0.1, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0), limits = c(-3,3)) +
        #geom_smooth(size=0.5) +
        xlab("Position") + ylab("Log R ratio") +
        facet_grid(Sample_ID~Chr, scales = "free", space = "free") +
        sams_pub_theme() +
        theme(panel.spacing.x = unit(0, "lines")) +
        theme(panel.border = element_rect(fill = NA))
        
    return(gg_log_ratio)
}
```

Plot BAF
```{r, cache=TRUE, warning=FALSE}
gg_baf <- plot_BAF(snp_fls = snp_fls, chrom_list = 1:22) +
    theme(axis.text.x = element_blank(), strip.text.y = element_text(angle = 0))
gg_baf
```

Plot Log Ratios
```{r, cache=TRUE, warning=FALSE}
gg_lr <- plot_LR(snp_fls = snp_fls, chrom_list = 1:22) +
    theme(axis.text.x = element_blank(), strip.text.y = element_text(angle = 0))
gg_lr
```

Write the plots to PNG files to make publication figures. There are too many data points for PDF files.
```{r, cache=TRUE}
png("snp_array/all_snp_array_baf_plots.png", width = 12, height = 8,
    units = "in", res = 300)
gg_baf
dev.off()

png("snp_array/all_snp_array_log_rato_plots.png", width = 12, height = 8,
    units = "in", res = 300)
gg_lr
dev.off()
```


For each sample, plot the chromosomes separately
```{r, cache=TRUE}

plot_LR_single_png <- function(x, chrom_list=chroms){
    
    snp_fl <- snp_fls[x]
    
    df <- read_snp(path = snp_fl, chroms = chrom_list)
    df$Chr <- factor(df$Chr, levels=c(1:22, "X", "Y"))

    gg_log_ratio <- ggplot(df, aes(x = Position, y = Log_R_Ratio)) +
        geom_point(size=0.1, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0), limits = c(-2,2)) +
        #geom_smooth(size=0.5) +
        xlab("Position") + ylab("Log R ratio") +
        facet_grid(Chr~., scales = "free", space = "free") +
        sams_pub_theme() +
        theme(panel.spacing.x = unit(0, "lines")) +
        theme(panel.border = element_rect(fill = NA))
    
    out <- str_c("snp_array/", names(snp_fls)[x], "_LR_contig_plots.png")
    
    ggsave(plot = gg_log_ratio, filename = out, device = "png",
           width = 12, height = 8, units = "in", dpi = 300, bg = "white")

}

lapply(X = 1:length(snp_fls), FUN = plot_LR_single_png)


```



```{r, cache=TRUE}
plot_BAFsingle_png <- function(x, chrom_list=chroms){
    
    snp_fl <- snp_fls[x]
    
    df <- lapply(snp_fl, read_snp, chroms=chrom_list) %>% do.call(rbind, .)
    df$Chr <- factor(df$Chr, levels=c(1:22, "X", "Y"))
    
    gg_freq <- ggplot(data = df, aes(y = B_Allele_Freq, x = Position)) +
        geom_point(size=0.1, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0)) +
        xlab("Position") + ylab("BAF") +
        facet_grid(Chr~., scales = "free", space = "free") +
        sams_pub_theme() +
        theme(panel.spacing.x = unit(0, "lines")) +
        theme(panel.border = element_rect(fill = NA))
    
    out <- str_c("snp_array/", names(snp_fls)[x], "_BAF_contig_plots.png")
    
    ggsave(plot = gg_freq, filename = out, device = "png",
           width = 12, height = 8, units = "in", dpi = 300, bg = "white")
}

lapply(X = 1:length(snp_fls), FUN = plot_BAFsingle_png)
```

Plot regions of interest
```{r, cache=TRUE}
plot_snp_region <- function(snp_fl, chrom, range){
    
    df <- lapply(snp_fl, read_snp, chroms=chrom) %>% do.call(rbind, .)
    df <- df[df$Position >= as.numeric(range[1]) & df$Position <= as.numeric(range[2]), ]
    
    gg_freq <- ggplot(data = df, aes(y = B_Allele_Freq, x = Position)) +
        geom_point(size=0.2, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0)) +
        xlab("Position") + ylab("BAF") +
        sams_pub_theme()
    gg_freq <- gg_freq + ggtitle(str_c(basename(snp_fl), "   chr", chrom, ":",
                                       as.integer(range[1]), "-",
                                       as.integer(range[2])))
    
    gg_log_ratio <- ggplot(df, aes(x = Position, y = Log_R_Ratio)) +
        geom_point(size=0.2, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0)) +
        geom_smooth(size=0.5) +
        xlab("Position") + ylab("Log R ratio") +
        sams_pub_theme()
    
    gg <- cowplot::plot_grid(plotlist = list(gg_freq,gg_log_ratio), nrow = 2)
    return(gg)
}

pdf("snp_array/region_plots_chr6.pdf")
lapply(snp_fls, plot_snp_region, chrom=6, range=c(1e8, 1.25e8))
dev.off()


pdf("snp_array/region_plots_chr4.pdf")
lapply(snp_fls, plot_snp_region, chrom=4, range=c(0, 5e7))
dev.off()


```


### Perform CNV detection

```{r, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("copynumber")
```


```{r}
library(copynumber)
```

Make matricies of BAF and LRR
```{r, cache=TRUE}
dat <- read_snp(path = snp_fls[1])

# measure must be one of "Log_R_Ratio" or "B_Allele_Freq"

make_mat <- function(file_list, measure="Log_R_Ratio"){
  
  dat_list <- lapply(file_list, read_snp)
  
  get_measure <- function(x){
    dat <- dat_list[[x]]
    val <- dat[ ,..measure]
    return(val)
    }
  
  dat_measure <- lapply(1:length(dat_list), get_measure)
  dat_measure <- do.call(cbind, dat_measure)
  colnames(dat_measure) <- names(file_list)
  
  loci <- data.frame(chr=dat_list[[1]]$Chr, pos=dat_list[[1]]$Position)
  
  df <- cbind(loci, dat_measure)
  
  df_sorted <- df[order(df$chr, df$pos, decreasing = FALSE), ]
  
  return(df_sorted)
}

mat_lrr <- make_mat(file_list = snp_fls, measure = "Log_R_Ratio")
mat_baf <- make_mat(file_list = snp_fls, measure = "B_Allele_Freq")

```


Filter rows with missing data
```{r, cache=TRUE}
keep <- complete.cases(mat_lrr)
mat_lrr <- mat_lrr[keep, ]
mat_baf <- mat_baf[keep, ]
```


Normalisation using  winsorize
```{r, message=FALSE, cache=TRUE}
lrr_win <- winsorize(mat_lrr, assembly = "hg19")
#baf_win <- winsorize(mat_baf, assembly = "hg19")
```

```{r, cache=TRUE}
plotGamma(mat_lrr, sample = 3, chrom=4, cex=3, cv = TRUE)
```

Call CNV segments
```{r, cache=TRUE}
allele_seg <- aspcf(logR = lrr_win, BAF = mat_baf, kmin = 10, assembly = "hg19")
```

Plot segmentation
```{r, cache=TRUE}
plotAllele(logR = lrr_win, BAF = mat_baf, segments = allele_seg, sample=c(1:6),chrom=c(20), ylim=c(-1,1))
```


```{r, cache=TRUE}

plotHeatmap(segments = allele_seg, upper.lim = 0.25)

pdf("snp_array/test.pdf", height = 16)
plotAberration(segments=allele_seg, thres.gain = 0.2,
               thres.loss = -0.25, layout = c(11,2), chrom = 1:22)
dev.off()
```




```{r, cache=TRUE}
allele_seg$width <- allele_seg$end.pos - allele_seg$start.pos
allele_seg[allele_seg$logR.mean < -0.25 & allele_seg$n.probes > 10, ]
```

Plot segment statistics
```{r, cache=TRUE}
gg_seg_size <- ggplot(data = allele_seg, aes(x = log(width))) +
  geom_density() +
  facet_grid(sampleID~.)
gg_seg_size
```

```{r, cache=TRUE}
gg_seg_lrr <- ggplot(data = allele_seg, aes(x = logR.mean)) +
  geom_density() +
  facet_grid(sampleID~.)
gg_seg_lrr


gg_seg_stat <- ggplot(data = allele_seg, aes(x = logR.mean, y = BAF.mean)) +
  geom_point() +
  facet_wrap(sampleID~.) +
  geom_vline(xintercept = c(-.2, .2)) +
  geom_hline(yintercept = 0.7) +
  xlim(c(-0.5, 0.5)) +
  ylim(c(0.5, 1))
gg_seg_stat
```

Plot InDel hits as matrix
```{r, cache=TRUE}
# Set heatmap thresholds
lrr_threshold <- 0.25
baf_threshold <- 0.7
#snp_count <- 10

allele_seg_filt <- allele_seg[(abs(allele_seg$logR.mean) > lrr_threshold) &
                                (allele_seg$BAF.mean > baf_threshold), ]

hits_mat <- matrix(data = 0, nrow = length(unique(allele_seg$sampleID)),
                   ncol = length(unique(allele_seg$chrom)))

rownames(hits_mat) <- c("Primed_32F_P33", "Primed_38F_P22",
                        "TNT_32F_P3", "TNT_38F_P3",
                        "NtP_32F_P17_P23", "NtP_38F_P12_P16")

colnames(hits_mat) <- c(1:22, "X")

hit_samples <- allele_seg_filt$sampleID %>% unique()


for (x in 1:length(hit_samples)) {
  
  hits_mat[hit_samples[x],
           allele_seg_filt$chrom[allele_seg_filt$sampleID == hit_samples[x]]] <- 1
}

pdf("snp_array/indel_hit_heatmap.pdf", width = 6, height = 1.5)
pheatmap(hits_mat, color = c("white", "red"),
         cluster_rows = FALSE, cluster_cols = FALSE, legend = FALSE, angle_col = 0)
dev.off()

pheatmap(hits_mat, color = c("white", "red"),
         cluster_rows = FALSE, cluster_cols = FALSE, legend = FALSE, angle_col = 0)
```

Plot indel segments
```{r, cache=TRUE}

logR = lrr_win, BAF = mat_baf, segments = allele_seg

plot_snp_region <- function(segment=allele_seg_filt[1, ],
                            lrr=lrr_win, baf=mat_baf, flank_multiple=2){
    
    plot_chrom <- segment$chrom
    segment_width <- segment$end.pos - segment$start.pos
    plot_range <- c((segment$start.pos - (segment_width*flank_multiple)),
                      segment$end.pos + (segment_width*flank_multiple))
      
    lrr_sub <- lrr[(lrr$chrom == plot_chrom) &
                     (lrr$pos > min(plot_range)) &
                     (lrr$pos < max(plot_range)), ]
    
    lrr_sub <- reshape2::melt(lrr_sub, id.vars=c("chrom", "pos"))
    colnames(lrr_sub)[1] <- "chr"
    lrr_sub$measure <- "LRR"
    
    baf_sub <- baf[(baf$chr == plot_chrom) &
                     (baf$pos > min(plot_range)) &
                     (baf$pos < max(plot_range)), ]
      
    baf_sub <- reshape2::melt(baf_sub, id.vars=c("chr", "pos"))
    baf_sub$measure <- "BAF"
    
    df <- rbind(lrr_sub, baf_sub)
    df$variable <- factor(df$variable,
                          levels=c("Primed_32F_P33", "Primed_38F_P22",
                        "TNT_32F_P3", "TNT_38F_P3",
                        "NtP_32F_P17_P23", "NtP_38F_P12_P16"))
    
    gg <- ggplot(data = df, aes(y = value, x = pos)) +
        geom_point(size=0.1, alpha=0.5) +
        scale_x_continuous(expand=c(0,0)) +
        scale_y_continuous(expand=c(0,0)) +
        facet_grid(measure~variable, scales = "free_y", space = "free") +
        xlab("Position") + 
        sams_pub_theme()
    
    return(gg)
}

pdf("snp_array/indel_hit_plots.pdf", width = 8, height = 6)
cowplot::plot_grid(
  plot_snp_region(segment = allele_seg_filt[1, ]),
  plot_snp_region(segment = allele_seg_filt[2, ]),
  plot_snp_region(segment = allele_seg_filt[3, ]), nrow = 3, align = "v")
dev.off()
```