4.2-CNV_analysis.Rmd

---
title: "CNV_analysis"
author: "Massimo Andreatta & Laura Yerly"
output: html_document
date: "2024-02-06"
---

```{r setup, include=FALSE,fig.width=16,fig.height=12}
renv::restore()
library(Seurat)
library(infercnv)
```

# Set path

```{r}
file_clean <- "cache/WHexp_annotated.rds"
seu <- readRDS(file_clean)
```

```{r}
table(seu$Sample, seu$patient_bcc)
table(seu$Sample, seu$annotation)
```

# Isolate tumor cells, normal keratinocytes and T cells for CNV analysis
```{r}
sub <- subset(seu, annotation %in% c("Tcell","Cancer_cells","Normal_Kerat"))
tab <- table(sub$Sample, sub$annotation)
tab
```
Focus on samples with enough cells
```{r eval=F}
ids <- tab[,"Cancer_cells"] >= 50 | tab[,"Normal_Kerat"] >= 50
pass <- names(ids)[ids]

sub <- subset(sub, subset=Sample %in% pass)
```

Downsample, for tests
```{r}
sub.list <- SplitObject(sub, split.by = "Sample")

sub.list.ds <- lapply(sub.list, function(x) {
  Idents(x) <- "annotation"
  subset(x, downsample = 2000)  #limit number of cells for largest group
})

sub.ds <- Reduce(merge, sub.list.ds)
table(sub.ds$Sample, sub.ds$annotation)
```


# InferCNV analysis

Generate gene position file from cellranger gtf: <https://github.com/broadinstitute/inferCNV/wiki/instructions-create-genome-position-file> got the python script on infercnv github and ran on Rstudio Terminal: python3 ../scripts/gtf_to_position_file.py --attribute_name gene_name /export/scratch/twyss/FKuonen_group/BCC_scRNAseq/refdata-gex-GRCh38-2020-A/genes/genes.gtf refdata-gex-GRCh38-2020-A_gen_pos.txt

```{r}
gene_order_file.path<-"_aux/refdata-gex-GRCh38-2020-A_gen_pos.txt"
outdir <- "cache/inferCNV_BCC_all_WH_TumKerT"
options("Seurat.object.assay.version" = "v3")

annot <- sub.ds@meta.data[,c("annotation","Sample")]

infercnv <- CreateInfercnvObject(raw_counts_matrix = sub.ds[["RNA"]]$counts,
                                 annotations_file = annot,
                                 gene_order_file = gene_order_file.path ,
                                 ref_group_names = c("Tcell"))
options(scipen = 100)
infercnv <- infercnv::run(infercnv,
                          cutoff=0.1, #threshold on gene expression
                          out_dir=outdir,
                          cluster_by_groups = FALSE,
                          denoise=TRUE,
                          leiden_resolution = 10^(-10),
                          hclust_method = "ward.D2",
                          num_threads = 8,
                          HMM=TRUE)
```


Interpret results
```{r eval=T}
scores=apply(infercnv@expr.data, 2, function(x){ sum(x < 0.9 | x > 1.1)/length(x) })

sub.ds$CNVpct <- scores

VlnPlot(sub.ds, features=c("CNVpct"), group.by = "annotation")
```

Run InferCNV separately for each sample
```{r}
seu.list <- SplitObject(sub, split.by = "patient_bcc")

infercnv.list <- lapply(names(seu.list), function(x) {
  
  obj <- seu.list[[x]]
  
  Idents(sub) <- "annotation"

  outdir <- sprintf("cache/inferCNV_%s_WH_TumKerT",x)
  annot <- as.data.frame(obj$annotation)
  
  infercnv.this <- CreateInfercnvObject(raw_counts_matrix = obj[["RNA"]]$counts,
                                   annotations_file = annot,
                                   gene_order_file = gene_order_file.path ,
                                   ref_group_names = c("Tcell"))
  
  options(scipen = 100)
  infercnv.this <- infercnv::run(infercnv.this,
                          cutoff=0.1, #threshold on gene expression
                          out_dir=outdir,
                          cluster_by_groups = FALSE,
                          denoise=TRUE,
                          leiden_resolution = 10^(-10),
                          hclust_method = "ward.D2",
                          num_threads = 8,
                          HMM=TRUE)
  infercnv.this
})
```