Single-Sample CNA Analysis

library(SCEVAN)
library(ggplot2)
library(dplyr)

# Load pre-processed MGH106 glioblastoma data
load(url("https://www.dropbox.com/s/b9udpvhnc2ez9pc/MGH106_data.RData?raw=1"))

# Examine data structure
cat("Dimensions:", dim(count_mtx), "\n")
cat("Genes:", nrow(count_mtx), "\n")
cat("Cells:", ncol(count_mtx), "\n")

# Preview
count_mtx[1:5, 1:5]

results <- pipelineCNA(
  count_mtx,
  sample = "MGH106",
  par_cores = 4,           # Adjust based on your system
  organism = "human",
  SUBCLONES = TRUE,        # Enable subclone detection
  beta_vega = 0.5,         # Segmentation parameter
  ClonalCN = TRUE,         # Infer clonal profile
  plotTree = TRUE          # Generate phylogenetic tree
)

# View classification results
head(results)

# Summary statistics
table(results$class)

# Subclone distribution (if detected)
if("subclone" %in% colnames(results)) {
  table(results$subclone, useNA = "ifany")
}

# List all output files
list.files("./output", pattern = "MGH106")

# Load the CNA matrix for custom analyses
load("./output/MGH106_CNAmtx.RData")
load("./output/MGH106_count_mtx_annot.RData")

dim(CNA_mtx_relat)

# Define a genomic region of interest (e.g., chromosome 7)
chr7_cn <- apply(
  CNA_mtx_relat[count_mtx_annot$seqnames == 7, ], 
  2, 
  mean
)

# Create data frame for plotting
cn_df <- data.frame(
  cell = names(chr7_cn),
  chr7_cn = chr7_cn,
  class = results[names(chr7_cn), "class"]
)

# Visualize
ggplot(cn_df, aes(x = class, y = chr7_cn, fill = class)) +
  geom_boxplot() +
  theme_minimal() +
  labs(
    title = "Chromosome 7 Copy Number by Cell Type",
    x = "Cell Classification",
    y = "Mean CN Ratio"
  ) +
  scale_fill_manual(values = c("tumor" = "#E74C3C", "normal" = "#3498DB"))

# Identify genes in amplified regions
chr7_genes <- rownames(count_mtx_annot)[count_mtx_annot$seqnames == 7]

# Compare expression between tumor subclones
if("subclone" %in% colnames(results)) {
  tumor_cells <- rownames(results)[results$class == "tumor"]
  
  # Get subclone assignments
  subclone_1 <- rownames(results)[results$subclone == 1 & !is.na(results$subclone)]
  subclone_2 <- rownames(results)[results$subclone == 2 & !is.na(results$subclone)]
  
  # Calculate fold changes for chr7 genes
  if(length(subclone_1) > 0 & length(subclone_2) > 0) {
    expr_s1 <- rowMeans(count_mtx[chr7_genes, subclone_1])
    expr_s2 <- rowMeans(count_mtx[chr7_genes, subclone_2])
    
    fc <- log2((expr_s1 + 1) / (expr_s2 + 1))
    
    # Top differentially expressed genes
    head(sort(fc, decreasing = TRUE), 10)
  }
}

# Specify known normal cell barcodes
known_normals <- c("cell_1", "cell_2", "cell_3")

results <- pipelineCNA(
  count_mtx,
  sample = "MGH106_custom",
  norm_cell = known_normals,
  FIXED_NORMAL_CELLS = FALSE  # Still allow automatic detection
)

# Coarser segmentation for noisy data
results_coarse <- pipelineCNA(
  count_mtx,
  sample = "MGH106_coarse",
  beta_vega = 1.0  # More conservative segmentation
)

# Define custom normal cell signatures
custom_signatures <- list(
  Astrocytes = c("GFAP", "AQP4", "SLC1A2", "SLC1A3"),
  Oligodendrocytes = c("MBP", "MOG", "PLP1", "OLIG2"),
  Microglia = c("CX3CR1", "P2RY12", "TMEM119", "AIF1")
)

results <- pipelineCNA(
  count_mtx,
  sample = "MGH106_custom_sig",
  AdditionalGeneSets = custom_signatures,
  SCEVANsignatures = TRUE  # Also use built-in signatures
)

# 1. Initial run with defaults
results_default <- pipelineCNA(count_mtx, sample = "sample_v1")

# 2. Review outputs and classifications
table(results_default$class)

# 3. Adjust parameters if needed
results_tuned <- pipelineCNA(
  count_mtx, 
  sample = "sample_v2",
  beta_vega = 0.7  # Adjusted based on review
)

# 4. Final analysis
results_final <- pipelineCNA(
  count_mtx,
  sample = "sample_final",
  SUBCLONES = TRUE,
  plotTree = TRUE
)