Skip to content

GitLab

Commit 79e28964 authored by Stefano Beretta's avatar Stefano Beretta
Browse files

Upload New File

parent 2fd8c621
Show whitespace changes
Inline Side-by-side
scripts/3_Visualization_Export.R 0 → 100644
View file @ 79e28964
library(SeuratObject)
library(Seurat)
library(harmony)
library(grid)
library(cowplot)
library(gridExtra)
library(RColorBrewer)
library(patchwork)
library(scales)
library(ggplot2)
library(openxlsx)
library(SingleR)
########################
### General Settings ###
########################
# Working dir
wdir <- "squadrito_livertumor2022_scrnaseq"
# Results dir
out_dir <- paste(wdir, "results", sep = "/")
dir.create(path = out_dir, showWarnings = F)
# Load Full object
Full_obj <- readRDS(paste(out_dir, "Full_DBR", "Full_DBR_final_DBR_labeled.rds", sep = "/"))
NKT_obj <- readRDS(paste(out_dir, "TNK", "TNK_DBR_final.rds", sep = "/"))
APC_obj <- readRDS(paste(out_dir, "APCs", "APCs_DBR_final.rds", sep = "/"))
# TCR recomputation
NKT_obj@meta.data$Samples <- "Undefined"
NKT_obj@meta.data$Samples[NKT_obj@meta.data$orig.ident=="GG-22"] <- "Combo_22"
NKT_obj@meta.data$Samples[NKT_obj@meta.data$orig.ident=="GG-11"] <- "Ctrl_11"
NKT_obj@meta.data$Samples[NKT_obj@meta.data$orig.ident=="GG-18"] <-"IFNa_18"
NKT_obj@meta.data$Samples[NKT_obj@meta.data$orig.ident=="GG-23"] <-"IFNa_23"
NKT_obj@meta.data$Samples[NKT_obj@meta.data$orig.ident=="GG-25"] <-"Ctrl_25"
NKT_obj@meta.data$Newlabels_detailed[NKT_obj@meta.data$Newlabels_detailed=="CD8 Teff_GG cells"] <-"CD8 Teff_ex cells"
NKT_obj@meta.data$Newlabels_detailed[NKT_obj@meta.data$Newlabels_detailed=="CD8 Teff_TK cells"] <-"CD8 Teff_ex cells"
###New ColFreq
NKT_obj@meta.data$IDS <- rownames(NKT_obj@meta.data)
df.VDJ <- NKT_obj@meta.data %>% select(IDS,Samples,cdr3_aa2) %>% filter(!is.na(cdr3_aa2)) %>%
group_by(Samples,cdr3_aa2) %>% mutate(ClonoFreq2=n())
NKT_obj@meta.data <- merge(NKT_obj@meta.data,df.VDJ,by.x=0,by.y=1,all.x=T)
rownames(NKT_obj@meta.data) <- NKT_obj@meta.data$IDS
###Make contingecy tables
#cloneType
NKT_obj@meta.data$cloneType <- "NA"
a1 <- 25
a2 <- 15
a3 <- 5
a4 <- 1
#cloneType
NKT_obj@meta.data$cloneType <- "NA"
# Assign values based on conditions
NKT_obj@meta.data$cloneType[NKT_obj@meta.data$ClonoFreq2 > a1 ] <- paste0("Hyperexpanded (x>",a1,")")
NKT_obj@meta.data$cloneType[NKT_obj@meta.data$ClonoFreq2 > a2 & NKT_obj@meta.data$ClonoFreq2 <= a1] <- paste0("Large (",a2,"< x <=",a1,")")
NKT_obj@meta.data$cloneType[NKT_obj@meta.data$ClonoFreq2 > a3 & NKT_obj@meta.data$ClonoFreq2 <= a2] <- paste0("Medium (",a3,"< x <=",a2,")")
NKT_obj@meta.data$cloneType[NKT_obj@meta.data$ClonoFreq2 > a4 & NKT_obj@meta.data$ClonoFreq2 <= a3] <- paste0("Small(",a4,"< x <=",a3,")")
NKT_obj@meta.data$cloneType[NKT_obj@meta.data$ClonoFreq2 == a4 ] <- paste0("Unique(x=",a4,")")
#####################################
## Visualisation and data analysis ##
#####################################
# Figure 8F
# Upper pannel
DimPlot(object = NKT_obj, reduction = "umap", group.by = "Newlabels_detailed",
pt.size = 1, label = T, split.by = "RNA_Group")+ylim(-6,7)
# Lower pannel
DimPlot(object = NKT_obj, reduction = "umap", group.by = "cloneType",
cols = c("brown","brown1","darkorange","grey","darkslategray","chartreuse4"),
pt.size = 1, label = F, split.by = "RNA_Group")+ylim(-6,7)
# Data analysis for Figure 8G (number of cells in cloneType)
df1 <- table(NKT_obj@meta.data$cloneType, NKT_obj@meta.data$Samples.x)
write.xlsx(df1,"TNK_TCR_clonotypes.xlsx")
# Data analysis for Figure 8H (number of cells in cloneType)
df1 <- table(NKT_obj@meta.data$ClonoFreq2,
NKT_obj@meta.data$Samples.x)/as.numeric(rownames(table(NKT_obj@meta.data$ClonoFreq2,
NKT_obj@meta.data$Samples.x)))
colSums(df1[as.numeric(rownames(df1))>1,])
# Figure S8E
DimPlot(object = APC_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1.2, label = T, label.size = 4, split.by = "RNA_Group") +
xlim(-6,6)
# Figure S8F
APC_obj@meta.data$RNA_Group <- factor(APC_obj@meta.data$RNA_Group, levels = rev(unique(APC_obj@meta.data$RNA_Group)))
MHCI_genes <- c("H2-T22", "H2-T23", "H2-D1","B2m", "Tap1", "Tap2", "Tapbp", "Psmb8", "Psmb9", "Cd86")
MHCII_genes <- c("H2-Aa","H2-Ab1", "H2-Eb1","H2-DMb1","H2-Oa", "Ciita", "Cd74", "Cd40")
IL10_gens <- c("Tgfb1", "Cebpb", "Il4ra", "Socs3", "Ccl24")
CTRL_genes <- c("Mmp8", "Tmem176b","Trem2", "Fn1")
DotPlot(object = APC_obj, features = c(MHCI_genes, MHCII_genes, IL10_gens, CTRL_genes),
group.by = "RNA_Group", dot.scale = 15, scale = F, dot.min = 0.1) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
# Figure S8G
# Feature plots
plist <- FeaturePlot(TNK_obj, reduction = "umap", features = c( "Pdcd1", "Havcr2","Ifng","Entpd1","Itga1","Cxcr6", "Prf1"),
pt.size = 1.2, ncol = 4, order = TRUE, combine=FALSE)
plist <- lapply(plist, function(g) {
g + ylim(c(-6,7))
})
CombinePlots(plist, ncol = 4)
# Dimplot
DimPlot(object = TNK_obj, reduction = "umap", group.by = "Newlabels_detailed", pt.size = 1.2, label = T, label.size = 4) + ylim(c(-6,7))
# Figure S8H
# Generate obj containing T cells with detected TCR
TNK_obj <- SetIdent(object = TNK_obj, value = "cloneType")
# Create DotPlot
TNK_VDJ_obj <- subset(TNK_obj, idents = c("NA"), invert = T)
DotPlot(TNK_VDJ_obj,
features = c("Entpd1", "Cd27", "Tnfrsf9", "Mir155hg", "Batf", "Tigit", "Pdcd1", "Havcr2", "Il7r",
"Ptprc", "Prf1", "Ifng", "Tnf", "Itga1", "Cd69"),
group.by = "cloneType", dot.scale = 11) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
####################
## Export results ##
####################
data_dir <- paste(wdir, "data", sep = "/")
dir.create(path = data_dir, showWarnings = F)
## Full ##
# Counts
full_counts <- Full_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "Full_DBR_final_DBR_labeled_counts.csv.gz", sep = "/"), "w")
write.csv(full_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
full_md <- Full_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "Full_DBR_final_DBR_labeled_metadata.csv.gz", sep = "/"), "w")
write.csv(x = full_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
full_umap <- Full_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "Full_DBR_final_DBR_labeled_umap.csv.gz", sep = "/"), "w")
write.csv(x = full_umap, gz_out_umap)
close(gz_out_umap)
## Markers Full object Newlabels
DefaultAssay(object = Full_obj) <- "RNA"
annotation <- "Newlabels"
Full_obj <- SetIdent(Full_obj, value = annotation)
tableFull_obj <- FindAllMarkers(Full_obj)
write.table(tableFull_obj,
paste(data_dir, "scRNAseq_DIFgenes_Full_obj_FC0.25 (Table S6).txt", sep = "/"),
quote = FALSE, sep = "\t")
## TNK ##
# Counts
TNK_counts <- TNK_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "TNK_DBR_final_counts.csv.gz", sep = "/"), "w")
write.csv(TNK_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
TNK_md <- TNK_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "TNK_DBR_final_metadata.csv.gz", sep = "/"), "w")
write.csv(x = TNK_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
TNK_umap <- TNK_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "TNK_DBR_final_umap.csv.gz", sep = "/"), "w")
write.csv(x = TNK_umap, gz_out_umap)
close(gz_out_umap)
## Markers TNK_obj Newlabels_detailed
DefaultAssay(object = TNK_obj) <- "RNA"
annotation <- "Newlabels_detailed"
TNK_obj <- SetIdent(TNK_obj, value = annotation)
tableTNK_obj <- FindAllMarkers(TNK_obj)
write.table(tableTNK_obj,
paste(data_dir, "scRNAseq_DIFgenes_TNK_FC0.25 (Table S8).txt", sep = "/"),
quote = FALSE, sep = "\t")
## APCs ##
# Counts
APC_counts <- APC_obj@assays$RNA@counts
gz_out_counts <- gzfile(paste(data_dir, "APCs_DBR_final_counts.csv.gz", sep = "/"), "w")
write.csv(APC_counts, gz_out_counts)
close(gz_out_counts)
# Meta Data
APC_md <- APC_obj@meta.data
gz_out_md <- gzfile(paste(data_dir, "APCs_DBR_final_metadata.csv.gz", sep = "/"), "w")
write.csv(x = APC_md, gz_out_md)
close(gz_out_md)
# UMAP Coordinates
APC_umap <- APC_obj@reductions$umap@cell.embeddings
gz_out_umap <- gzfile(paste(data_dir, "APCs_DBR_final_umap.csv.gz", sep = "/"), "w")
write.csv(x = APC_umap, gz_out_umap)
close(gz_out_umap)
## Markers APC_obj Newlabels_detailed
DefaultAssay(object = APC_obj) <- "RNA"
annotation <- "Newlabels_detailed"
APC_obj <- SetIdent(APC_obj, value = annotation)
tableAPC_obj <- FindAllMarkers(APC_obj)
write.table(tableAPC_obj,
paste(data_dir, "scRNAseq_DIFgenes_APC_FC0.25 (Table S7).txt", sep = "/"),
quote = FALSE, sep = "\t")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment