CNV Calling

Introduction

Copy Number Variation (CNV) analysis can be performed using various tools that use multiple approaches to define CNV profile from WES or WGS data. Usually CNV analysis is done on two steps, the first is CNV calling where the CNV profile of each locus of the genomic data is defined, then the second step involves segmentation of the CNV calls to define segments of amplification and deletion events and potential break points. Recommended readings on different approaches for CNV calling and segmentation:

• Copy Number Variation (CNV) Analysis: Complete Guide
• Copy Number Variation (CNV) Calling
• Whole-genome CNV analysis: advances in computational approaches (Pirooznia et al. 2015)
• On the core segmentation algorithms of copy number variation detection tools (Zhang et al. 2024)

We will use CNVKit(Talevich et al. 2016), which uses read depth as a proxy for CNV estimation, which is the best method used for Illumina short read sequencing. CNVKit is a python package that can run through bash. It has extensive documentation here for usage and different plotting functionalities.

CNVKit requires the essential resources for the CNV calling as follows:

• Reference genome fasta file
• Target regions bed file : regions captured by the WES kit (It’s recommended to use the one provided by the exome capture kit manufacturer)
• RefFlat gene annotation file

We have downloaded these files in the Data_Preprocessing step and placed them in the data/references folder as follows

ImportantDownloading References

Human (hg38) references

### Load necessary modules
ml bwa
ml samtools
ml gatk
ml picard-tools

### Create reference directory
mkdir data/references
refDir=data/references

### Download References Genome from UCSC (Ignore these lines if you already downloaded the references genome)
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz -O .${refDir}/hg38.fa.gz
gunzip ${refDir}/hg38.fa.gz

### Index References Genome
bwa index ${refDir}/hg38.fa
samtools faidx ${refDir}/hg38.fa
gatk CreateSequenceDictionary -R ${refDir}/hg38.fa

### Download Resources from GATK Google Cloud Bucket (Recommneded for GATK workflow)
## https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0;tab=objects?prefix=&forceOnObjectsSortingFiltering=false
## https://console.cloud.google.com/storage/browser/gatk-best-practices/somatic-hg38;tab=objects?prefix=&forceOnObjectsSortingFiltering=false
## NOTE : You might need to download crcmod via `pip install crcmod`
gsutil cp gs://genomics-public-data/resources/broad/hg38/v0/Homo_sapiens_assembly38.* $refDir

wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/refFlat.txt.gz -O ${refDir}/refFlat.txt.gz
gunzip ${refDir}/refFlat.txt.gz

wget https://www.twistbioscience.com/sites/default/files/resources/2022-12/hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.bed -O $refDir/hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.bed

Mouse (mm39) references

### Load necessary modules
ml bwa
ml samtools
ml gatk
ml picard-tools

### Create reference directory
mkdir data/references
refDir=data/references

########################################
### Download References Genome from UCSC
########################################
wget https://hgdownload.soe.ucsc.edu/goldenPath/mm39/bigZips/mm39.fa.gz -O ${refDir}/mm39.fa.gz
gunzip ${refDir}/mm39.fa.gz

###########################
### Index References Genome
###########################
bwa index ${refDir}/mm39.fa
samtools faidx ${refDir}/mm39.fa
gatk CreateSequenceDictionary -R ${refDir}/mm39.fa

###########################
### refFlat mm39
###########################
wget https://hgdownload.soe.ucsc.edu/goldenPath/mm39/database/refFlat.txt.gz -o ${refDir}/refFlat.txt.gz
gunzip ${refDir}/refFlat.txt.gz

###########################
### WES Capture bed file
###########################
## WES capture kit bed file is downloaded as follows
## (https://kb.10xgenomics.com/hc/en-us/articles/115004150923-Where-can-I-find-the-Agilent-Target-BED-files-)
## Since the bait file is based on mm9, liftover tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver) was used to
## map the coordinates to mm39

wget --timestamping \
    'ftp://hgdownload.soe.ucsc.edu/goldenPath/mm9/liftOver/mm9ToMm39.over.chain.gz' \
    -O ${refDir}/mm9ToMm39.over.chain.gz

gunzip ${refDir}/mm9ToMm39.over.chain.gz

tail -n +3 ${refDir}/S0276129/S0276129_Regions.bed > ${refDir}/S0276129/S0276129_Regions_no_header.bed

./liftOver \
${refDir}/S0276129/S0276129_Regions_no_header.bed \
${refDir}/mm9ToMm39.over.chain \
${refDir}/S0276129/S0276129_Regions_mm39.bed unMapped

CNV Analysis Pipeline

1. CNVKit Analysis

• Batchmode
• Step-by-step analysis
• Genemetrics

2. CNV Calling Results Visualization

• Whole Genome CNV Profile Heatmap
• Gene-specific CNV Heatmaps
• Top Amplified and Deleted Genes Heatmap
• Specific genes log2ratio and copy number value from genemetrics output
• Calculating gene-level log2 ratio using weighted bin-level log2 ratio
• CNV mutation burden
• Plot Specific Gene Coverage and CNV using karyoploteR

CNVKit Analysis

Batchmode

CNVKit can be run in batchmode, where all samples are processed in one command. This is useful when you have multiple tumor and normal samples. It uses all normal samples to create a pooled reference rather than using the paired Tumor-Normal analysis. CNVKit authors recommend using this approach

To analyze a cohort sequenced on a single platform, we recommend combining all normal samples into a pooled reference, even if matched tumor-normal pairs were sequenced – our benchmarking showed that a pooled reference performed slightly better than constructing a separate reference for each matched tumor-normal pair. Furthermore, even matched normals from a cohort sequenced together can exhibit distinctly different copy number biases (see Plagnol et al. 2012 and Backenroth et al. 2014); reusing a pooled reference across the cohort provides some consistency to help diagnose such issues.

TipAdding ‘chr’ prefix to references

Sometimes the reference genome fasta file and the bed files have different chromosome naming conventions. For example, the reference genome might have chromosome names like chr1, chr2, etc., while the bed files might have chromosome names like 1, 2, etc. This mismatch can lead to errors during CNVKit analysis. To resolve this, you can use sed or awk to remove the ‘chr’ prefix as needed as follows:

ref=data/references/Homo_sapiens_assembly38.fasta
bait=data/references/hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.bed
refFlat=data/references/refFlat.txt

sed 's/^>chr/>/' $ref > data/references/Homo_sapiens_assembly38_wout_chr.fasta
sed 's/^chr//' $bait > data/references/hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.wout_chr.bed
awk -F'\t' '{gsub(/^chr/, "", $3); print}' $refFlat > data/references/refFlat_w_chr.txt

We also showed how to add ‘chr’ to samples bam files in Data_Preprocessing

We use the following command to run CNVKit

cnvkit.py access <ref_genome> -o <ref_genome_access_file>

cnvkit.py batch <tumor_samples> \
--normal <normal_samples> \
--targets <capture_bed_file> \
--annotate <refFlat_file> \
--fasta <ref_genome> \
--access <ref_genome_access_file> \
--output-reference <ref_genome_cnn_output> \
--output-dir <output_dir> \
--diagram --scatter

You can also use CNVKit without control/normal samples as follows. Review Biostars post and GitHub issue

cnvkit.py batch <tumor_samples> \
    -normal \
    -targets <capture_bed_file> \
    -fasta <ref_genome> \
    --access <ref_genome_access_file> \
    --output-reference <ref_genome_cnn_output>

Note that batch runs multiple steps at once, you can review the detailed steps and run them one-by-one for more customization as demonstrated in the CNVKit documentation

NoteSLURM Script

It’s recommended to run CNVKit on HPC cluster when you have large number of samples. Here is a slurm script to run CNVKit in batchmode. You might want to install CNVKit in a conda environment as follows: conda create -n cnvkit -c bioconda cnvkit then activate it using conda activate cnvkit. This script requires a manifest file CNVKit_manifest.txt that contains three columns, the first is the sample name, the second is the bam file path, the third is the sample type (Tumor/Normal). The results will be saved in the outputs/CNVKit directory.

column1	column2	column3
sample1	path/to/sample1.bam	Tumor
sample2	path/to/sample2.bam	Normal
…	…	…

#!/bin/bash

#SBATCH --account=
#SBATCH --job-name='CNVKit_batchmode'
#SBATCH --output=logs/CNVKit_batchmode.log
#SBATCH --partition=standard
#SBATCH --mem=128G
#SBATCH --cpus-per-task=32
#SBATCH --time=24:00:00

source activate ~/.conda/envs/cnvkit/
outputDir=outputs/CNVKit
mkdir -p ${outputDir}
ref=data/references/Homo_sapiens_assembly38.fasta ## or Homo_sapiens_assembly38_w_chr.fasta
access=${outputDir}/access.hg38.bed
bait=data/references/hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.bed ## or hg38_exome_v2.0.2_targets_sorted_validated.re_annotated.w_chr.bed
refFlat=data/references/refFlat.txt ## or data/references/refFlat_w_chr.txt
samples_manifest=CNVKit_manifest.txt
normal_samples=$(cat $samples_manifest | grep Normal | cut -f9 | tr '\n' ' ')
tumor_samples=$(cat $samples_manifest | grep -v Normal | cut -f9 | tr '\n' ' ')

cnvkit.py access $ref -o $access

cnvkit.py batch $tumor_samples \
--normal $normal_samples \ ## Set to nothing if no normal samples
--targets $bait \
--annotate $refFlat \
--fasta $ref \
--access $access \
--output-reference $outputDir/hg38_reference.cnn \
--output-dir $outputDir \
--diagram --scatter

TipStep-by-step analysis

To run CNVKit step-by-step rather than using batch

#| echo: false

## If you would like to run it on sample-wise approach, where each tumor sample is compared to the paired normal sample, use the following script

#!/bin/bash

#SBATCH --account=
#SBATCH --job-name='CNVKit_analysis'
#SBATCH --output=logs/CNVKit_analysis_%a.log

#SBATCH --partition=standard
#SBATCH --mem=16G
#SBATCH --cpus-per-task=8
#SBATCH --time=08:00:00

#SBATCH --array=1-<number_of_samples>

source activate ~/.conda/envs/cnvkit/
bam_dir=outputs/GATK_preprocessing
outputDir=outputs/CNVKit
mkdir -p ${outputDir}
target=${outputDir}/hg38.target.bed
antitarget=${outputDir}/hg38.antiarget.bed
samples_manifest=CNVKit_manifest.txt
normal_reference=${outputDir}/normal_reference.cnn
tumor_sample=$(cat $samples_manifest | grep -v 'Normal' | sed -n ${SLURM_ARRAY_TASK_ID}p | cut -f1)
bam_path=$(cat $samples_manifest | grep -v 'Normal' | sed -n ${SLURM_ARRAY_TASK_ID}p | cut -f8)
sample_bam=$(echo ${bam_dir}/$(basename $bam_path .bam)_sorted.bam)

cnvkit.py coverage ${sample_bam} ${target} \
-o ${outputDir}/${tumor_sample}.targetcoverage.cnn

cnvkit.py coverage ${sample_bam} ${antitarget} \
-o ${outputDir}/${tumor_sample}.antitargetcoverage.cnn

cnvkit.py fix ${outputDir}/${tumor_sample}.targetcoverage.cnn \
${outputDir}/${tumor_sample}.antitargetcoverage.cnn \
${normal_reference} \
-o ${outputDir}/${tumor_sample}.cnr

cnvkit.py segment ${outputDir}/${tumor_sample}.cnr \
-o ${outputDir}/${tumor_sample}.cns

cnvkit.py call ${outputDir}/${tumor_sample}.cns \
-o ${outputDir}/${tumor_sample}.call.cns

cnvkit.py export seg ${outputDir}/${tumor_sample}.call.cns -o ${outputDir}/${tumor_sample}.seg

cnvkit.py bintest -a 1 -o ${outputDir}/${tumor_sample}_bintest.cnr \
${outputDir}/${tumor_sample}.cnr

cnvkit.py call ${outputDir}/${tumor_sample}_genemetrics.cnr \
-o ${outputDir}/${tumor_sample}_genemetrics_calls.cnr

CNVKit Outputs

CNVKit provides several outputs that you can use for downstream analysis.

• .cnr provided bin-level log2 ratio, that represents a metric for amplification (if positive) or deletion (if negative) at each bin (where the bins are defined during the CNVKit analysis)
• .cns provided segmented log2 ratio, where the bins are merged into segments of amplification or deletion events. This should be more reliable than the bin-level data.
• .call.cns provided segmented log2 ratio with discrete copy number calls (e.g., gain, loss, neutral)
• Other outputs can be found in the CNVKit documentation here

Genemetrics

CNVKit provides a genemetrics command that summarizes CNV calls at the gene level. This can be useful for identifying genes that are frequently amplified or deleted across samples. You can use call command after genemetrics to classify the CNV calls into discrete categories (e.g., gain, loss, neutral).

cnvkit.py genemetrics <sample_cnr_output> \
-t 0 \
-s <sample_cns_output> \
-o <sample_genemetrics_output>

cnvkit.py call <sample_genemetrics_output> \
-o <sample_genemetrics_output_calls>

To run this for all sample, we can use the following SLURM script that needs outputDir and samples_manifest variables to be defined as in the previous scripts.

NoteSLURM Script

#!/bin/bash

#SBATCH --account=
#SBATCH --job-name='CNVKit_batchmode'
#SBATCH --output=logs/CNVKit_batchmode.log
#SBATCH --partition=standard
#SBATCH --mem=128G
#SBATCH --cpus-per-task=32
#SBATCH --time=24:00:00
#SBATCH --array=1-<number_of_samples>

outputDir=outputs/CNVKit
samples_manifest=CNVKit_manifest.txt
tumor_sample=$(cat $samples_manifest | grep -v 'Normal' | sed -n ${SLURM_ARRAY_TASK_ID}p | cut -f1)

cnvkit.py genemetrics ${outputDir}/${tumor_sample}.cnr \
-t 0 \
-s ${outputDir}/${tumor_sample}.cns \
-o ${outputDir}/${tumor_sample}_genemetrics.cnr

cnvkit.py call ${outputDir}/${tumor_sample}_genemetrics.cnr \
-o ${outputDir}/${tumor_sample}_genemetrics_calls.cnr

CNV Calling Results Visualization

Whole Genome CNV Profile Heatmap

cnvkit.py heatmap -d -o outputs/CNVKit/cnv_heatmap_cns.png outputs/CNVKit/*.call.cns

Gene-specific CNV Heatmaps

You have to specify the genomic coordinates of the gene of interest using -c parameter. Here are examples for Pten and Trp53 genes in mouse genome (mm39)

## Pten 
cnvkit.py heatmap -d -o outputs/CNVKit/pten_heatmap.png outputs/CNVKit/*.cnr -c chr19:32734897-32803560

## Trp53 
cnvkit.py heatmap -d -o outputs/CNVKit/Trp53_heatmap.png outputs/CNVKit/N*.cnr -c  chr11:69471185-69482698

Top Amplified and Deleted Genes Heatmap

library(tidyverse)
library(pheatmap)

genemetrics <- list.files(
  "outputs/CNVKit/",
  pattern = "genemetrics_calls.cnr",
  full.names = T
)
genemetrics <- lapply(genemetrics, read.table, header = T)
names(genemetrics) <- gsub(
  "genemetrics_calls.cnr",
  "",
  list.files("outputs/CNVKit/", pattern = "genemetrics_calls.cnr")
)
genemetrics <- bind_rows(genemetrics, .id = 'sample_id') %>%
  mutate(event = ifelse(cn > 2, "gain", ifelse(cn < 2, "loss", "Neutral")))

freq <- table(genemetrics$gene, genemetrics$event) %>%
  data.frame

top_amp_genes <- freq %>%
  filter(Var2 == 'gain' & Freq != 0) %>%
  filter(!grepl(x = .$Var1, pattern = ".*Rik$")) %>%
  arrange(desc(Freq)) %>%
  distinct(Var1) %>%
  pull(Var1) %>%
  as.character()

top_del_genes <- freq %>%
  filter(Var2 == 'loss' & Freq != 0) %>%
  filter(!grepl(x = .$Var1, pattern = ".*Rik$")) %>%
  arrange(desc(Freq)) %>%
  distinct(Var1) %>%
  pull(Var1) %>%
  as.character()

Specific genes log2ratio and copy number value from genemetrics output

gois <- c("Trp53", "Pten")
plot_df <- genemetrics %>% filter(gene %in% gois) %>% complete(sample_id, gene)
plot_df$log2[is.na(plot_df$log2)] <- 0
log2 <- ggplot(plot_df, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  geom_text(aes(label = round(p_bintest, 2))) +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"))

plot_df$cn[is.na(plot_df$cn)] <- 2
cn <- ggplot(plot_df, aes(sample_id, gene)) +
  geom_tile(aes(fill = factor(cn)), colour = "white") +
  geom_text(aes(label = round(p_bintest, 2))) +
  scale_fill_manual(values = c('2' = 'white', '1' = 'blue', '0' = 'darkblue')) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"))

cowplot::plot_grid(log2, cn, nrow = 2)

Calculating gene-level log2 ratio using weighted bin-level log2 ratio

suppressPackageStartupMessages({
  library(pheatmap)
  library(janitor)
  library(tidyverse)
  library(readxl)
  library(GenomicRanges)
})

samples_info <- read_xlsx(
  "../data/Mouse Sample Key_11 13 2023 .xlsx",
  sheet = 'For Ahmed'
) %>%
  as.data.frame() %>%
  clean_names() %>%
  filter(wes_rna_or_both != 'RNA') %>%
  mutate(
    spon_vs_cell_line_vs_cldt = ifelse(
      tumor_vs_control %in% c("Normal liver", "Normal fat"),
      tumor_vs_control,
      spon_vs_cell_line_vs_cldt
    )
  ) %>%
  dplyr::select(c(
    'sample_id',
    "mouse_id",
    "sex",
    "tumor_vs_control",
    "spon_vs_cell_line_vs_cldt"
  ))

tumor_samples_names <- samples_info %>%
  filter(!tumor_vs_control %in% c("Normal_fat", "Normal_liver", "AAV8_fat")) %>%
  pull(sample_id)
spont_tumor <- samples_info %>%
  filter(
    !tumor_vs_control %in% c("Normal_fat", "Normal_liver", "AAV8_fat") &
      spon_vs_cell_line_vs_cldt == 'Spont Tumor p53PTEN'
  ) %>%
  arrange(desc(tumor_vs_control)) %>%
  pull(sample_id)
cellLine <- samples_info %>%
  filter(
    !tumor_vs_control %in% c("Normal_fat", "Normal_liver", "AAV8_fat") &
      spon_vs_cell_line_vs_cldt != 'Spont Tumor p53PTEN'
  ) %>%
  arrange(desc(tumor_vs_control)) %>%
  pull(sample_id)

tumor_cnr <- lapply(
  list.files(
    "outputs/CNVKit_revised_run/",
    pattern = '[0-9].cnr$',
    full.names = T
  ),
  read.delim,
  header = T,
  sep = "\t"
) %>%
  `names<-`(gsub(
    ".cnr",
    "",
    list.files("outputs/CNVKit_revised_run/", pattern = '[0-9].cnr$')
  ))

pten_exon_5_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  pten_gr <- GRanges("chr19:32777261-32777499")
  pten_overlap <- findOverlaps(cnr, pten_gr)
  pten_cnr <- cnr[queryHits(pten_overlap), ]
  log2 <- weighted.mean(pten_cnr$log2, pten_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Pten_Exon5') %>%
  dplyr::select(sample_id, gene, log2, event)

# pten_logR <- lapply(tumor_samples_names, function(x){
#   cnr <- GRanges(tumor_cnr[[x]])
#   pten_gr <- GRanges("chr19:32734897-32803560")
#   pten_overlap <- findOverlaps(cnr, pten_gr)
#   pten_cnr <- cnr[queryHits(pten_overlap),]
#   log2 <- weighted.mean(pten_cnr$log2, pten_cnr$weight)
#   event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
#   return(data.frame(log2, event))
# }) %>% `names<-`(tumor_samples_names) %>%
#   bind_rows(.id = 'sample_id') %>%
#   mutate(gene = 'Pten') %>%
#   dplyr::select(sample_id, gene, log2, event)

trp53_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  trp53_gr <- GRanges("chr11:69471185-69482699")
  trp53_overlap <- findOverlaps(cnr, trp53_gr)
  trp53_cnr <- cnr[queryHits(trp53_overlap), ]
  log2 <- weighted.mean(trp53_cnr$log2, trp53_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Trp53') %>%
  dplyr::select(sample_id, gene, log2, event)

plot_df <- rbind(pten_exon_5_logR, trp53_logR)

plot_df_spont <- plot_df %>%
  filter(sample_id %in% spont_tumor) %>%
  mutate(sample_id = factor(sample_id, levels = spont_tumor))

log2_spont <- ggplot(plot_df_spont, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

plot_df_celline <- plot_df %>%
  filter(sample_id %in% cellLine) %>%
  mutate(sample_id = factor(sample_id, levels = cellLine))
log2_cellline <- ggplot(plot_df_celline, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

mdm2_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  mdm2_gr <- GRanges("chr10:117524780-117546663")
  mdm2_overlap <- findOverlaps(cnr, mdm2_gr)
  mdm2_cnr <- cnr[queryHits(mdm2_overlap), ]
  log2 <- weighted.mean(mdm2_cnr$log2, mdm2_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Mdm2') %>%
  dplyr::select(sample_id, gene, log2, event)
frs2_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  frs2_gr <- GRanges("chr10:116905332-116984415")
  frs2_overlap <- findOverlaps(cnr, frs2_gr)
  frs2_cnr <- cnr[queryHits(frs2_overlap), ]
  log2 <- weighted.mean(frs2_cnr$log2, frs2_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Frs2') %>%
  dplyr::select(sample_id, gene, log2, event)

hmga2_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  hmga2_gr <- GRanges("chr10:120197180-120312374")
  hmga2_overlap <- findOverlaps(cnr, hmga2_gr)
  hmga2_cnr <- cnr[queryHits(hmga2_overlap), ]
  log2 <- weighted.mean(hmga2_cnr$log2, hmga2_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Hmga2') %>%
  dplyr::select(sample_id, gene, log2, event)

cdk4_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(tumor_cnr[[x]])
  cdk4_gr <- GRanges("chr10:126899403-126903789")
  cdk4_overlap <- findOverlaps(cnr, cdk4_gr)
  cdk4_cnr <- cnr[queryHits(cdk4_overlap), ]
  log2 <- weighted.mean(cdk4_cnr$log2, cdk4_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Cdk4') %>%
  dplyr::select(sample_id, gene, log2, event)

plot_df <- rbind(mdm2_logR, frs2_logR, hmga2_logR, cdk4_logR)

plot_df_spont <- plot_df %>%
  filter(sample_id %in% spont_tumor) %>%
  mutate(sample_id = factor(sample_id, levels = spont_tumor))
log2_spont <- ggplot(plot_df_spont, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

plot_df_celline <- plot_df %>%
  filter(sample_id %in% cellLine) %>%
  mutate(sample_id = factor(sample_id, levels = cellLine))
log2_cellline <- ggplot(plot_df_celline, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    axis.title.y = element_blank(),
    axis.title.x = element_blank()
  )

# event <- ggplot(plot_df, aes(sample_id, gene)) +
#   geom_tile(aes(fill = factor(event)), colour = "white") +
#   scale_fill_manual(values = c('neutral' = 'white', 'gain' = 'red', 'loss' = 'darkblue')) +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"))

::::::::::::::::::::::


library(pheatmap)
library(janitor)
library(tidyverse)
library(readxl)

samples_info <- read_xlsx(
  "../data/Mouse Sample Key_11 13 2023 .xlsx",
  sheet = 'For Ahmed'
) %>%
  as.data.frame() %>%
  clean_names() %>%
  filter(wes_rna_or_both != 'RNA') %>%
  mutate(
    spon_vs_cell_line_vs_cldt = ifelse(
      tumor_vs_control %in% c("Normal liver", "Normal fat"),
      tumor_vs_control,
      spon_vs_cell_line_vs_cldt
    )
  ) %>%
  dplyr::select(c(
    'sample_id',
    "mouse_id",
    "sex",
    "tumor_vs_control",
    "spon_vs_cell_line_vs_cldt"
  ))

tumor_samples_names <- samples_info %>%
  filter(!tumor_vs_control %in% c("Normal fat", "Normal liver")) %>%
  pull(sample_id)

genemetrics <- list.files(
  "outputs/CNVKit/",
  pattern = "_genemetrics_calls.cnr",
  full.names = T
)
genemetrics <- lapply(genemetrics, read.table, header = T)
names(genemetrics) <- gsub(
  "_genemetrics_calls.cnr",
  "",
  list.files("outputs/CNVKit/", pattern = "_genemetrics_calls.cnr")
)
genemetrics <- bind_rows(genemetrics, .id = 'sample_id') %>%
  mutate(event = ifelse(cn > 2, "gain", ifelse(cn < 2, "loss", "Neutral")))

freq <- table(genemetrics$gene, genemetrics$event) %>%
  data.frame
# reshape(idvar = 'Var1', timevar = 'Var2', direction = "wide") %>%
# `colnames<-`(c("gene","gain","loss","neutral"))

top_amp_genes <- freq %>%
  filter(Var2 == 'gain' & Freq != 0) %>%
  filter(!grepl(x = .$Var1, pattern = ".*Rik$")) %>%
  arrange(desc(Freq)) %>%
  distinct(Var1) %>%
  pull(Var1) %>%
  as.character()

top_del_genes <- freq %>%
  filter(Var2 == 'loss' & Freq != 0) %>%
  filter(!grepl(x = .$Var1, pattern = ".*Rik$")) %>%
  arrange(desc(Freq)) %>%
  distinct(Var1) %>%
  pull(Var1) %>%
  as.character()


pten_exon_5_logR <- lapply(tumor_samples_names, function(x) {
  cnr <- GRanges(read.delim(paste0("outputs/CNVKit/", x, ".cnr")))
  pten_gr <- GRanges("chr19:32777261-32777499")
  pten_overlap <- findOverlaps(cnr, pten_gr)
  pten_cnr <- cnr[queryHits(pten_overlap), ]
  log2 <- weighted.mean(pten_cnr$log2, pten_cnr$weight)
  event <- ifelse(log2 < -0.25, "loss", ifelse(log2 > +0.2, "gain", "neutral"))
  return(data.frame(log2, event))
}) %>%
  `names<-`(tumor_samples_names) %>%
  bind_rows(.id = 'sample_id') %>%
  mutate(gene = 'Pten_Exon5') %>%
  select(sample_id, gene, log2, event)

gois <- c("Trp53", "Pten")
plot_df <- genemetrics %>%
  filter(gene %in% gois) %>%
  complete(sample_id, gene) %>%
  select(sample_id, gene, log2, event)
plot_df$log2[is.na(plot_df$log2)] <- 0
plot_df$event[is.na(plot_df$event)] <- "neutral"
plot_df <- rbind(plot_df, pten_exon_5_logR)

log2 <- ggplot(plot_df, aes(sample_id, gene)) +
  geom_tile(aes(fill = log2), colour = "white") +
  scale_fill_gradient2(
    low = "blue",
    mid = 'white',
    high = "red",
    midpoint = 0
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"))

event <- ggplot(plot_df, aes(sample_id, gene)) +
  geom_tile(aes(fill = factor(event)), colour = "white") +
  scale_fill_manual(
    values = c('neutral' = 'white', 'gain' = 'red', 'loss' = 'darkblue')
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"))

cowplot::plot_grid(log2, event, nrow = 2)

CNV mutation burden

setwd(dirname(rstudioapi::getSourceEditorContext()$path))

library(tidyverse)
library(janitor)
library(GenomicRanges)


samples_info <- readxl::read_xlsx(
  "../data/Mouse Sample Key_11 13 2023 .xlsx",
  sheet = 'For Ahmed'
) %>%
  clean_names()

cns_files <- list.files(
  "outputs/CNVKit_revised_run/",
  pattern = 'call.cns$',
  full.names = T
)
cns <- lapply(cns_files, read.table, header = T) %>%
  `names<-`(gsub(
    ".call.cns",
    "",
    list.files("outputs/CNVKit_revised_run/", pattern = 'call.cns$')
  ))

samples_info <- samples_info %>% filter(sample_id %in% names(cns))
samples_info <- samples_info %>% filter(tumor_vs_control != "Normal fat")

spont_tumor <- samples_info %>%
  filter(spon_vs_cell_line_vs_cldt == 'Spont Tumor p53PTEN') %>%
  pull(sample_id)

cell_lines <- samples_info %>%
  filter(spon_vs_cell_line_vs_cldt != 'Spont Tumor p53PTEN') %>%
  pull(sample_id)

cns_spont_tumor <- cns[spont_tumor]

cns_spont_tumor_mutational_burden <- lapply(cns_spont_tumor, function(x) {
  freq <- data.frame(table(x$cn)) %>%
    filter(Var1 != 2) %>%
    `colnames<-`(c("CN", "Freq")) %>%
    mutate(CN = as.numeric(as.character(CN))) %>%
    mutate(Event = ifelse(CN < 2, "Loss", "Gain"))
}) %>%
  bind_rows(.id = "sample")

ggplot(
  cns_spont_tumor_mutational_burden,
  aes(x = sample, y = Freq, fill = Event)
) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggtitle("No. of Unbalanced Segments (Spontaneous Tumor)")

cns_spont_cellLine <- cns[cell_lines]

cns_cellLine_mutational_burden <- lapply(cns_spont_cellLine, function(x) {
  freq <- data.frame(table(x$cn)) %>%
    filter(Var1 != 2) %>%
    `colnames<-`(c("CN", "Freq")) %>%
    mutate(CN = as.numeric(as.character(CN))) %>%
    mutate(Event = ifelse(CN < 2, "Loss", "Gain"))
}) %>%
  bind_rows(.id = "sample")

ggplot(
  cns_cellLine_mutational_burden,
  aes(x = sample, y = Freq, fill = Event)
) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ggtitle("No. of Unbalanced Segments (Cell Line)")


## Calculating percent altered
cns_fil <- lapply(cns, function(x) {
  x %>% filter(cn != 2)
})
cns_gr <- lapply(cns_fil, GRanges)
cns_gr <- lapply(cns_gr, function(x) {
  x$fraction <- (width(x) / 2.7e9)
  return(x)
})
percent_altered <- lapply(cns_gr, function(x) {
  sum(x$fraction)
}) %>%
  unlist() %>%
  data.frame() %>%
  rownames_to_column() %>%
  `colnames<-`(c("sample_id", "percent_altered"))
metadata <- samples_info %>%
  filter(spon_vs_cell_line_vs_cldt == 'Spont Tumor p53PTEN') %>%
  select(sample_id, tumor_vs_control, ti_ls)
percent_altered <- merge(percent_altered, metadata, by = 'sample_id')
percent_altered$condition <- paste0(
  percent_altered$tumor_vs_control,
  " (",
  percent_altered$ti_ls,
  ")"
)

ggplot(percent_altered, aes(x = condition, y = percent_altered)) +
  geom_boxplot() +
  geom_point() +
  xlab("") +
  ylab("% Genome Altered") +
  theme_bw()

percent_altered <- percent_altered[order(percent_altered$percent_altered), ]
percent_altered$sample_id <- factor(
  percent_altered$sample_id,
  levels = percent_altered$sample_id
)
ggplot(
  percent_altered,
  aes(x = sample_id, y = percent_altered, fill = condition)
) +
  geom_bar(stat = 'identity') +
  xlab("") +
  ylab("% Genome Altered") +
  labs(fill = "Grade") +
  coord_flip()

percent_altered$condition <- factor(
  percent_altered$condition,
  levels = c("WDLPS (Low TILs)", "DDLPS (High TILs)", "DDLPS (Low TILs)")
)
percent_altered <- percent_altered[order(percent_altered$condition), ]
percent_altered$sample_id <- factor(
  percent_altered$sample_id,
  levels = percent_altered$sample_id
)
ggplot(
  percent_altered,
  aes(x = sample_id, y = percent_altered, fill = condition)
) +
  geom_bar(stat = 'identity') +
  xlab("") +
  ylab("% Genome Altered") +
  labs(fill = "Grade") +
  coord_flip()

Plot Specific Gene Coverage and CNV using karyoploteR

Generate gene bed

First we will need to generate bed files for the coverage of the genes of interest from the bam file using bedtools. We will need the same manifest files that has a column for sample ID and bam file path. We also need to specify the genes of interest coordinates, for example here we are interested in Trp53 and Pten genes in mouse genome (mm39) and set this using trp53_coord="chr11:69468307-69485577" & pten_coord="chr19:32717731-32820726"

#!/bin/bash

#SBATCH --account=
#SBATCH --job-name='generate_gene_bed_$a'
#SBATCH --output=logs/generate_gene_bed_%a.log
#SBATCH --partition=standard
#SBATCH --mem=2G
#SBATCH --cpus-per-task=1
#SBATCH --time=00:15:00
#SBATCH --array=1-<number_of_samples>

ml bedtools2

samples_manifest=outputs/samples_manifest.txt
sample=$(cat $samples_manifest | cut -f1 | uniq | sed -n ${SLURM_ARRAY_TASK_ID}p)
bam_file_path=$(cat $samples_manifest | cut -f2 | uniq | sed -n ${SLURM_ARRAY_TASK_ID}p)
outputDir=outputs/bed_files
mkdir -p $outputDir

trp53_coord="chr11:69468307-69485577"
samtools view -b $bam_file_path $trp53_coord > ${sample}_trp53_temp.bam
pten_coord="chr19:32717731-32820726"
samtools view -b $bam_file_path $pten_coord > ${sample}_pten_temp.bam

bedtools bamtobed -i ${sample}_trp53_temp.bam > ${outputDir}/${sample}.bed
bedtools bamtobed -i ${sample}_pten_temp.bam >> ${outputDir}/${sample}.bed

rm ${sample}_trp53_temp.bam ${sample}_pten_temp.bam

Plot Gene-specific Coverage and CNV Profile

Here we will need to load the coverage bed files and the CNVKit cnr files into R and use karyoploteR to plot the coverage and CNV profile for the genes of interest. We will use Trp53 gene as an example here. We will need a vector for sample names that can be imported from the samples_manifest.txt file used above.

library(karyoploteR)
library(regioneR)
library(zoo)
library(janitor)
library(readxl)
library(TxDb.Mmusculus.UCSC.mm39.knownGene)
library(TxDb.Hsapiens.UCSC.hg39.knownGene)
library(tidyverse)
library(fields)

samples_info <- read.table('samples_manifest.txt', header = F, sep = "\t") %>%
  as.data.frame()

coverage <- lapply(samples_info[,1], function(x) {
  message(x)
  gr <- GRanges(
    read.delim(
      paste0("outputs/bed_files/", x, ".bed"),
      sep = "\t",
      header = F
    ) %>%
      dplyr::select(V1, V2, V3) %>%
      `colnames<-`(c("chromosome", "start", "end"))
  )
}) %>%
  `names<-`(samples_info[,1])

# Importing CNVKit cnr files
cnr <- lapply(samples_info[,1],
  function(x) {
    gr <- GRanges(read.delim(
      paste0("outputs/CNVKit/", x, ".cnr"),
      sep = "\t",
      header = T
    ))
  }) %>%
  `names<-`(samples_info[,1])

target_bed <- GRanges(
  read.table("data/references/S0276129_Regions_mm39.target.bed") %>%
    `colnames<-`(c("chromosome", "start", "end", "gene"))
)

cnr <- lapply(cnr, function(x) {
  GenomicRanges::merge(x, target_bed)
})

all_cnr <- bind_rows(lapply(cnr, data.frame), .id = 'sample_id')

tumor_samples <- c() ## Here you have to specify the tumor samples names vector

tumor_cnr <- all_cnr %>% filter(sample_id %in% tumor_samples) %>% GRanges()
tumor_trp53_overlaps <- findOverlaps(
  spont_tumor_cnr,
  GRanges("chr11:69468307-69485577")
)
tumor_trp53_cnr <- tumor_cnr[
  queryHits(tumor_trp53_overlaps),
]
tumor_ymin = min(tumor_trp53_cnr$log2)
tumor_ymax = max(tumor_trp53_cnr$log2)

Trp53 <- plotKaryotype(
  genome = 'mm39',
  plot.type = 2,
  chromosomes = c("chr11"),
  zoom = GRanges("chr11:69468307-69485577")
)
kpAddBaseNumbers(
  Trp53,
  tick.dist = 2000,
  tick.len = 10,
  add.units = T,
  cex = 0.5,
  minor.ticks = T,
  minor.tick.len = 5,
  minor.tick.dist = 100
)
Trp53_genes.data <- makeGenesDataFromTxDb(
  txdb = TxDb.Mmusculus.UCSC.mm39.knownGene,
  karyoplot = Trp53
)
Trp53_genes.data <- addGeneNames(Trp53_genes.data)
Trp53_genes.data <- mergeTranscripts(Trp53_genes.data)

total.tracks <- 42
kpPlotGenes(
  Trp53,
  Trp53_genes.data,
  gene.name.position = "left",
  data.panel = 2,
  r0 = 0.3,
  r1 = 0.7
)

at <- autotrack(current.track = 1, total.tracks = total.tracks, margin = 0.2)
kpRect(
  Trp53,
  chr = 'chr11',
  x0 = 69468307,
  x1 = 69485577,
  y0 = 1,
  y1 = 1.1,
  r0 = at$r0,
  r1 = at$r1,
  data.panel = 1,
  col = 'darkgrey',
  border = 'darkgrey'
)
at <- autotrack(current.track = 2, total.tracks = total.tracks, margin = 0.2)
kpPlotCoverage(
  Trp53,
  data = coverage[['N1307']],
  r0 = at$r0,
  r1 = at$r1,
  data.panel = 1
)
kpAddLabels(
  Trp53,
  data.panel = 1,
  r0 = at$r0,
  r1 = at$r1,
  labels = c("N1307"),
  cex = 0.5
)
at <- autotrack(current.track = 3, total.tracks = total.tracks, margin = 0.2)
sample_cnr_data <- tumor_trp53_cnr[
  spont_tumor_trp53_cnr$sample_id == "SAMPLE_NAME" ## You have to specify the sample name here
]
kpHeatmap(
  Trp53,
  data = sample_cnr_data,
  r0 = at$r0,
  r1 = at$r1,
  y = sample_cnr_data$log2,
  data.panel = 1,
  colors = color_palette,
  ymin = spont_tumor_ymin,
  ymax = spont_tumor_ymax
)

image.plot(
  legend.only = TRUE,
  zlim = c(spont_tumor_ymin, spont_tumor_ymax),
  col = color_palette,
  legend.args = list(text = "log2 ratio", side = 4, line = 2.5, cex = 1),
  smallplot = c(0.89, 0.9, 0.9, 0.98)
)

References

Pirooznia, Mehdi, Fernando S. Goes, and Peter P. Zandi. 2015. “Whole-Genome CNV Analysis: Advances in Computational Approaches.” Frontiers in Genetics Volume 6 - 2015. https://doi.org/10.3389/fgene.2015.00138.

Talevich, Eric, A Hunter Shain, Thomas Botton, and Boris C Bastian. 2016. “CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing.” PLoS Computational Biology 12 (4): e1004873.

Zhang, Yibo, Wenyu Liu, and Junbo Duan. 2024. “On the Core Segmentation Algorithms of Copy Number Variation Detection Tools.” Briefings in Bioinformatics 25 (2): bbae022. https://doi.org/10.1093/bib/bbae022.