Axi-cel.Rmd

---
title: "Patients-tisa-cell-10x"
output: html_document
date: "2023-02-28"
---

## load packages
```{r}
#load packages
package_list <- c(
  "tidyverse", "Seurat", "patchwork", "SeuratDisk", "ggplot2", "ggrepel", "multimode", "dplyr", "EnhancedVolcano", "writexl", "cowplot", "RColorBrewer", "sctransform", "ggpubr", "ggVennDiagram", "rstatix", "Signac", "Libra"
)
invisible(lapply(package_list, require, character.only = TRUE))

set.seed(35) # set a seed to ensure the analysis is reproducible 

# create a directory for the results
results_dir <- "C:/Single Cell Sequencing/3 Axi-cel patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/results"
if (!dir.exists(results_dir)) {
  dir.create(results_dir, recursive = TRUE)
}

```

## load 10x data
```{r}
P01 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac01/")
P02 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac02/")
P03 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac03/")
P04 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac04/")
P05 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac05/")
P06 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac06/")
P07 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac07/")
P08 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac08/")
P09 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac09/")
P10 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac10/")
P11 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac11/")
P12 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac12/")
P13 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac13/")
P14 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac14/")
P15 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac15/")
P16 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac16/")
P17 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac17/")
P18 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac18/")
P19 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac19/")
P20 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac20/")
P21 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac21/")
P22 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac22/")
P23 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac23/")
P24 <- Read10X(data.dir = "C:/Single Cell Sequencing/patient CAR T cell infusion products/R-analysis/patient CAR T cell infusion products/scRNA.data/ac24/")
```


##create seurat object 

```{r}
#Filter genes detected in fewer than three cells
P01 <- CreateSeuratObject(counts = P01, project = "P01", min.cells = 3)
P02 <- CreateSeuratObject(counts = P02, project = "P02", min.cells = 3)
P03 <- CreateSeuratObject(counts = P03, project = "P03", min.cells = 3)
P04 <- CreateSeuratObject(counts = P04, project = "P04", min.cells = 3)
P05 <- CreateSeuratObject(counts = P05, project = "P05", min.cells = 3)
P06 <- CreateSeuratObject(counts = P06, project = "P06", min.cells = 3)
P07 <- CreateSeuratObject(counts = P07, project = "P07", min.cells = 3)
P08 <- CreateSeuratObject(counts = P08, project = "P08", min.cells = 3)
P09 <- CreateSeuratObject(counts = P09, project = "P09", min.cells = 3)
P10 <- CreateSeuratObject(counts = P10, project = "P10", min.cells = 3)
P11 <- CreateSeuratObject(counts = P11, project = "P11", min.cells = 3)
P12 <- CreateSeuratObject(counts = P12, project = "P12", min.cells = 3)
P13 <- CreateSeuratObject(counts = P13, project = "P13", min.cells = 3)
P14 <- CreateSeuratObject(counts = P14, project = "P14", min.cells = 3)
P15 <- CreateSeuratObject(counts = P15, project = "P15", min.cells = 3)
P16 <- CreateSeuratObject(counts = P16, project = "P16", min.cells = 3)
P17 <- CreateSeuratObject(counts = P17, project = "P17", min.cells = 3)
P18 <- CreateSeuratObject(counts = P18, project = "P18", min.cells = 3)
P19 <- CreateSeuratObject(counts = P19, project = "P19", min.cells = 3)
P20 <- CreateSeuratObject(counts = P20, project = "P20", min.cells = 3)
P21 <- CreateSeuratObject(counts = P21, project = "P21", min.cells = 3)
P22 <- CreateSeuratObject(counts = P22, project = "P22", min.cells = 3)
P23 <- CreateSeuratObject(counts = P23, project = "P23", min.cells = 3)
P24 <- CreateSeuratObject(counts = P24, project = "P24", min.cells = 3)


```

## merge all seurat
```{r}
merged_10x <- merge(P01,
                     y =c(P02,P03,P04,P05,P06,P07,P08,P09,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,P23,P24), add.cell.ids = c("lib1","lib2","lib3","lib4","lib5","lib6","lib7","lib8","lib9","lib10","lib11","lib12","lib13","lib14","lib15","lib16","lib17","lib18","lib19","lib20","lib21","lib22","lib23","lib24"),
  project = "patient CART"
)
```

## QC and filtering

```{r}
merged_10x[["percent.mt"]] <- PercentageFeatureSet(merged_10x, pattern = "^MT-")
# this filtering criteria is based on the original publication for this dataset!
merged_10x_filtered <- subset(merged_10x, subset = nFeature_RNA > 200 & nFeature_RNA < 7000 & percent.mt < 15)

```

## Normalization

```{r}
CART.patients <- NormalizeData(merged_10x_filtered)
```

## data Frame

```{r}
data_df_10xCAR <- data.frame(Expression=CART.patients[["RNA"]]@data["FMC63-CD19SCFV",]) # CAR-10x
data_df_10xCD4 <- data.frame(Expression=CART.patients[["RNA"]]@data["CD4",]) # CD4-10x
data_df_10xCD8A <- data.frame(Expression=CART.patients[["RNA"]]@data["CD8A",]) # CD8A-10x

nmodes(data = data_df_10xCAR$Expression, bw = bw.nrd0(data_df_10xCAR$Expression), full.result = T) # CAR-10x
nmodes(data = data_df_10xCD4$Expression, bw = bw.nrd0(data_df_10xCD4$Expression), full.result = T) # CD4-10x
nmodes(data = data_df_10xCD8A$Expression, bw = bw.nrd0(data_df_10xCD8A$Expression), full.result = T) # CD8A-10x

#CAR
locmodes(data_df_10xCAR$Expression, mod0=7,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),display=T)
### Define which are the antimodes and add the number in the expm1() to calculate the molecules per cell to be used in log1p(mols/cell)
expm1(0.1530759) 

VlnPlot(CART.patients, features = "FMC63-CD19SCFV", pt.size = 0.01, group.by = "orig.ident")+geom_hline(yintercept = log1p(0.1654134))+NoLegend()

#CD4
locmodes(data_df_10xCD4$Expression, mod0=6,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),display=T)
### Define which are the antimodes and add the number in the expm1() to calculate the molecules per cell to be used in log1p(mols/cell)
expm1(0.1198674)

VlnPlot(CART.patients, features = "CD4", pt.size = 0.01, group.by = "orig.ident")+geom_hline(yintercept = log1p(0.1273474))+NoLegend()

#CD8

locmodes(data_df_10xCD8A$Expression, mod0=2,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),display=T)
### Define which are the antimodes and add the number in the expm1() to calculate the molecules per cell to be used in log1p(mols/cell)
expm1(0.1787789) 

VlnPlot(CART.patients, features = "CD8A", pt.size = 0.01, group.by = "orig.ident")+geom_hline(yintercept = log1p(0.1957563))+NoLegend()

```


## CAR positive
```{r}
#CARrna

CAR <- WhichCells(object = CART.patients, expression = `FMC63-CD19SCFV` > log1p(0.1654134))

# cell labels are returned

head(CAR)

# we then add a new column to the meta data based on the presence of the cell label from WhichCells

CART.patients$CAR<- ifelse(colnames(CART.patients) %in% CAR, "CAR+", "CAR-")

# set the identities to the new column(s)
Idents(CART.patients) <- "CAR"


# CD4 RNA
CD4 <- WhichCells(object = CART.patients, expression = `CD4` > log1p(0.1273474))

# cell labels are returned

head(CD4)

# we then add a new column to the meta data based on the presence of the cell label from WhichCells

CART.patients$CD4<- ifelse(colnames(CART.patients) %in% CD4, "CD4+", "CD4-")

# set the identities to the new column(s)
Idents(CART.patients) <- "CD4"


# CD8 RNA
CD8 <- WhichCells(object = CART.patients, expression = `CD8A` > log1p(0.1957563))

# cell labels are returned

head(CD8)

# we then add a new column to the meta data based on the presence of the cell label from WhichCells

CART.patients$CD8<- ifelse(colnames(CART.patients) %in% CD8, "CD8+", "CD8-")

# set the identities to the new column(s)
Idents(CART.patients) <- "CD8"



# combine columns


CART.patients$CD4_CAR <- paste0(CART.patients$CD4 ,"_", CART.patients$CAR)

CART.patients$CD8_CAR <- paste0(CART.patients$CD8 ,"_", CART.patients$CAR)


CART.patients@meta.data


```

## Findmarkers
```{r}

#### CD4/CD8 CAR

Idents(CART.patients) <- "CD8_CAR"
CD8CAR <- FindMarkers(CART.patients, ident.1= "CD8+_CAR-", ident.2 = "CD8+_CAR+" , logfc.threshold = 0.25, min.pct = 0.2, test.use = "wilcox", assay = "RNA")
CD8CAR$FDR <- p.adjust(CD8CAR$p_val, method = "BH")
CD8CAR$gene = row.names(CD8CAR)
write_xlsx(CD8CAR, "results/DEG/CD4CAR.CD8CAR.xlsx", col_names = TRUE)
EnhancedVolcano(CD8CAR, lab = rownames(CD8CAR), x ="avg_log2FC", y ="FDR", title = "patients. CD8_CAR+ versus CD8_CAR-", FCcutoff = 0.25, labSize = 5.0, pointSize = 4, legendPosition = 'none', drawConnectors = T, max.overlaps = 30, pCutoff = 0.05, colAlpha = 1, cutoffLineType = 'twodash', cutoffLineWidth = 0.8, gridlines.major = F, gridlines.minor = F,labFace = "bold", xlim = c(-0.8,0.8), legendLabels = c("n.s", "|log2FC|>0.25", "FDR<0.05", "|log2FC|>0.25 & FDR<0.05"), selectLab = c("IFITM1", "IFITM2", "HLA-DQA2", "HLA-DQA1", "HLA-B"))


Idents(CART.patients) <- "CD4_CAR"
CD4CAR <- FindMarkers(CART.patients, ident.1= "CD4+_CAR-", ident.2 = "CD4+_CAR+" , logfc.threshold = 0.25, min.pct = 0.2, test.use = "wilcox", assay = "RNA")
CD4CAR$FDR <- p.adjust(CD4CAR$p_val, method = "BH")
CD4CAR$gene = row.names(CD4CAR)
write_xlsx(list(`CD8CAR` = CD8CAR, `CD4CAR` =CD4CAR), "results/DEG/CD4CAR.CD8CAR.xlsx", col_names = TRUE)
EnhancedVolcano(CD4CAR, lab = rownames(CD4CAR), x ="avg_log2FC", y ="FDR", title = "patients. CD4_CAR+ versus CD4_CAR-", FCcutoff = 0.25, labSize = 5.0,  labFace = "bold",pointSize = 4, legendPosition = 'none', drawConnectors = T, max.overlaps = 30, pCutoff = 0.05, colAlpha = 1, cutoffLineType = 'twodash', cutoffLineWidth = 0.8, gridlines.major = F, gridlines.minor = F,xlim = c(-0.8,0.8), legendLabels = c("n.s", "|log2FC|>0.25", "FDR<0.05", "|log2FC|>0.25 & FDR<0.05"), selectLab = c("IFITM1", "IFITM2", "HLA-A", "HLA-E", "HLA-B", "HLA-C", "ISG20"))

ggsave("results/vln/CD4CARLab5.pdf", width = 5, height = 5)

```

## violin plots

```{r}

# choose a vector of genes
gene <- c("IFITM1","IFITM2","IFITM3","STAT1","ISG20", "SAMHD1","IFI6", "ADAR", "OAS2", "MX1", "APOBEC3G", "DUSP2", "DUSP4", "CTLA4", "LAG3", "HAVCR2")


Idents(CART.patients) <- "patient_CART"
compare <- list( c("P01_CAR-", "P01_CAR+"),
                 c("P02_CAR-", "P02_CAR+"),
                 c("P03_CAR-", "P03_CAR+"),
                 c("P04_CAR-", "P04_CAR+"),
                 c("P05_CAR-", "P05_CAR+"),
                 c("P06_CAR-", "P06_CAR+"),
                 c("P07_CAR-", "P07_CAR+"),
                 c("P08_CAR-", "P08_CAR+"),
                 c("P09_CAR-", "P09_CAR+"),
                 c("P10_CAR-", "P10_CAR+"),
                 c("P11_CAR-", "P11_CAR+"),
                 c("P12_CAR-", "P12_CAR+"),
                 c("P13_CAR-", "P13_CAR+"),
                 c("P14_CAR-", "P14_CAR+"),
                 c("P15_CAR-", "P15_CAR+"),
                 c("P16_CAR-", "P16_CAR+"),
                 c("P17_CAR-", "P17_CAR+"),
                 c("P18_CAR-", "P18_CAR+"),
                 c("P19_CAR-", "P19_CAR+"),
                 c("P20_CAR-", "P20_CAR+"),
                 c("P21_CAR-", "P21_CAR+"),
                 c("P22_CAR-", "P22_CAR+"),
                 c("P23_CAR-", "P23_CAR+"),
                 c("P24_CAR-", "P24_CAR+"))




vln_df <- list()
plot_list_vln_df <- list()
pwc <- list()
p.val.size <- 6
plot.title.size <- 20
set.seed(42)
# change number according to total groups plotted
mycolors <- c( "darkred", "darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred","darkblue", "darkred", "darkblue")

for(i in gene){

#  violin plot without noise
vln_df[[i]] = data.frame(Expression=CART.patients[["RNA"]]@data[i,], Subset = CART.patients$patient_CAR
                         
                         )

# Statistics
pwc[[i]] <- vln_df[[i]] %>% 
  wilcox_test(Expression ~ Subset, p.adjust.method = "none", comparisons = compare)
pwc[[i]] <- pwc[[i]] %>% add_xy_position(x = "Subset")
pwc[[i]]$FDR <- p.adjust(pwc[[i]]$p, method = "BH")#, n=length(gene))
pwc[[i]] <- add_significance(pwc[[i]], p.col = "FDR", output.col = "FDR.signif")

# Graphs
# Add noise for optimal visualization purposes like VlnPlot() 
noise <- rnorm(n = length(x = vln_df[[i]][, "Expression"])) / 100000
vln_df[[i]]$Expression <- vln_df[[i]]$Expression  + noise

# Variable stats position.
### OBS! MIGHT NEED TO BE ADJUSTED for the hight of the p-value symbols and lines!!###
y.position <- seq(max(vln_df[[i]]$Expression),max(vln_df[[i]]$Expression)+ 0.5 , by=0.5)

# Variable ylim range
ylim <- ylim(0, max(y.position)+0.1)

# violin plot with noise
if(i != gene[length(gene)]){
  if (any(i == gene[seq(1, length(gene),by=7)])){
  plot_list_vln_df[[i]] <- vln_df[[i]] %>%
ggplot(mapping=aes(x=Subset, y=Expression)) + geom_violin(mapping=aes(x=Subset, y=Expression, fill=Subset),scale = "width")+scale_fill_manual(values = alpha((mycolors),0.9))+ geom_boxplot(outlier.size = 0.01, width=0.3, size=0.5, alpha=0.4) + 
    theme_bw()+
    xlab("") +
  ggtitle(label=i) +
  stat_pvalue_manual(pwc[[i]],  hide.ns = F, y.position = y.position, label = "FDR.signif", size = p.val.size, tip.length = 0) +
  labs(subtitle = get_test_label(vln_df[[i]] %>% kruskal_test(Expression ~ Subset) %>% mutate(KW_FDR = p.adjust (p, method='BH')) %>% add_significance(p.col = "KW_FDR"), p.col = "KW_FDR.signif", detailed=F, description = NULL))+ylim +
  theme(panel.grid = element_blank(),axis.title.y = element_text(size = 12), axis.text.y = element_text(size=12), axis.text.x = element_blank(), axis.ticks.x = element_blank(), plot.title = element_text(hjust = 0.5, size=plot.title.size, face = "bold"), plot.subtitle = element_blank())+ylab("Expression")+ theme(legend.position = "none" , aspect.ratio = 0.3)
} 
  else {
  plot_list_vln_df[[i]] <- vln_df[[i]] %>%
ggplot(mapping=aes(x=Subset, y=Expression)) + geom_violin(mapping=aes(x=Subset, y=Expression, fill=Subset),scale = "width")+scale_fill_manual(values = alpha((mycolors),0.9))+ geom_boxplot(outlier.size = 0.01, width=0.3, size=0.5, alpha=0.4) + 
    theme_bw()+
  xlab("") +
  ggtitle(label=i) +
  stat_pvalue_manual(pwc[[i]],  hide.ns = F, y.position =  y.position, label = "FDR.signif", size = p.val.size, tip.length = 0) +
  labs(subtitle = get_test_label(vln_df[[i]] %>% kruskal_test(Expression ~ Subset) %>% mutate(KW_FDR = p.adjust (p, method='BH')) %>% add_significance(p.col = "KW_FDR"), p.col = "KW_FDR.signif", detailed=F, description = NULL))+ylim +
  theme(panel.grid = element_blank(),axis.title.y = element_text(size = 12), axis.text.y = element_text(size=12), axis.text.x = element_blank(), axis.ticks.x = element_blank(), plot.title = element_text(hjust = 0.5, size=plot.title.size, face = "bold"), plot.subtitle = element_blank())+ylab("Expression")+ theme(legend.position = "none",  aspect.ratio = 0.3)
}
}
else{plot_list_vln_df[[i]] <- vln_df[[i]] %>%
ggplot(mapping=aes(x=Subset, y=Expression)) + geom_violin(mapping=aes(x=Subset, y=Expression, fill=Subset),scale = "width")+scale_fill_manual(values = alpha((mycolors),0.9))+ geom_boxplot(outlier.size = 0.01, width=0.3, size=0.5, alpha=0.4) + 
    theme_bw()+
  xlab("") +
  ggtitle(label=i) +
  stat_pvalue_manual(pwc[[i]],  hide.ns = F, y.position =  y.position, label = "FDR.signif", size = p.val.size, tip.length = 0) +
  labs(subtitle = get_test_label(vln_df[[i]] %>% kruskal_test(Expression ~ Subset) %>% mutate(KW_FDR = p.adjust (p, method='BH')) %>% add_significance(p.col = "KW_FDR"), p.col = "KW_FDR.signif", detailed=F, description = NULL))+ylim +
  theme(panel.grid = element_blank(),axis.title.y = element_text(size = 12), axis.text.y = element_text(size=12), axis.text.x = element_blank(), axis.ticks.x = element_blank(), plot.title = element_text(hjust = 0.5, size=plot.title.size, face = "bold"), plot.subtitle = element_blank()) +
  theme(legend.title = element_text(size = 12), legend.text = element_text(size = 12), legend.key.size = unit(0.3, "cm"), aspect.ratio = 0.3)+ylab("Expression")
}
}
```

```{r}
# Set the factor "f" to adjust the height of the exported pdf and play around with width
f =5
w = 15
col = 1
h = f*if(length(gene)==any(seq(col,length(gene),by=col))){length(gene)} else{(length(gene)/col)+0.5}
n = if(length(gene)==any(seq(col,length(gene),by=col))){length(gene)/col} else{(length(gene)/col)+0.5}

ggsave("results/vln/allpatients5.pdf", plot=ggarrange(plotlist=plot_list_vln_df, ncol = col, nrow=n, common.legend = T, legend = "top"), width = w, height = h, limitsize = F)
```