Skip to content

GitLab

Commit a78a5561 authored by Marco Monti's avatar Marco Monti
Browse files

I added the manuscript's scripts

parents
Expand all Show whitespace changes
Inline Side-by-side
CB2025_figure_1_RNAseq.R 0 → 100644
View file @ a78a5561
library("readxl")
library("dplyr")
#####Directories#####
us <- "/Users/Squadrito/"
us <-"C:/Users/bresesti.chiara/"
wdir1509<-paste0(us,"/Dropbox (HSR Global)/SquadritoM_1509_RNASeq_QIAseq_UPX/QIAseqUltraplexRNA_181342/primary_analysis")
wdir1510<-paste0(us,"/Dropbox (HSR Global)/SquadritoM_1510_RNA_miRNA_QIAseq_UPX")
fdir <- paste0(us,"/Dropbox (HSR Global)/CancerGeneTherapy/Cancer Gene Therapy/MS/2024 Bresesti et al/Scripts/plots and tables used in figures")
#### Figure 1B ####
setwd(wdir1509)
df1 <- as.data.frame(read_excel("QIAseqUltraplexRNA_181342.xlsx",sheet=3,col_names = T,skip = 1))#sheet: umis.genes.polyA-mouse
df1 <-df1[,-c(1,3:6)] #keep only gene name and UMI counts
df1[,-1] <-apply(df1[,-1],2,function(x){x/sum(x)*1000000}) #UMI normalized by CPM
gene_vector <- c("Cd19", "Ms4a1", "Fcer2a", "Ighm", "Cd8a", "Xcr1", "Itgae", "Itgax",
"Ccr2", "Itgam", "Mgl2", "Cd68", "Vcam1", "Csf1r", "Adgre1",
"Siglec1", "Hmox1", "Timd4", "Vsig4", "Clec4f", "Marco", "Pecam1",
"Tek", "Lyve1", "Stab2")
df2 <- df1[df1$gene%in%gene_vector,] %>% arrange(factor(gene, levels=gene_vector))
rownames(df2) <- df2[,1]
df2.scaled <- as.data.frame(t(scale(t(df2[-1])))) #Zscore normalization. Scaled only works on columns, so need to transform
setwd(fdir)
openxlsx::write.xlsx(df2.scaled,"Figure_1B_table.xlsx", rowNames=T)
#the df was exported and used to create a heatmap using Graphpad
#### Figure 1C ####
# "edgeR_DGE_res_volcano.pdf" (in the 'wdir1509' folder) was imported in illustrator
# and merged with the plot generated by the code below
setwd(wdir1510)
df1 <- as.data.frame(read_excel("173308.all_samples.summary.xlsx",sheet=2,col_names = T))#sheet: miRNA_piRNA
df1 <-df1[,1:7]#UMI
rownames(df1)<-df1$miRNA
df2 <- apply(df1[,-1], 2, function(x) log(x)) #Log counts
df2 <- as.data.frame(limma::normalizeCyclicLoess(df2, weights = NULL, span=0.7, iterations = 5, method = "pairs")) #Cyclic loess normalization
colnames(df2)<-c("B cell","RPM","cDC1","cDC2","LSEC","KC") #Rename columns
upper.panel<-function(x, y, ...){
points(((x+y)/2),(x-y), cex=0.6, pch=19, col="grey")
above <- (x-y-1)*((x+y)/2-5) > 5 & ((x+y)/2)>5
points(((x+y)/2)[above],(x-y)[above], col="red", cex=0.6, pch=19)
below <- (x-y+1)*((x+y)/2-5) < -5 & ((x+y)/2)>5
points(((x+y)/2)[below],(x-y)[below], col="blue", cex=0.6, pch=19)
} #function for MA plot
pairs(df2[,1:6], lower.panel = NULL, upper.panel = upper.panel,
ylim=c(-8.5,8.5), xlim=c(2,14), cex.labels = 2) #pairwise plot
#### Figure 1D ####
setwd(wdir1510)
df1 <- as.data.frame(read_excel("173308.all_samples.summary.xlsx",sheet=2,col_names = T))#sheet: miRNA_piRNA
df1 <- df1[-which(grepl("piR",df1[,1])),1:7]#UMI and piRNA removing
df1[,1] <- gsub("/.*","",df1[,1]) #leave only first miRNA for ambiguous entries
df1[,-1] <- apply(df1[,-1],2,function(x){x/sum(x,na.rm=T)*1000000})#UMI normalized by CPM
df.families <- as.data.frame(read_excel("Analysis CB/miRNA Family.xlsx",sheet=1,col_names = T))#miRNA families from miRBase
df.families <- df.families[df.families[,3]==10090,c(4,1)] #select mouse entries
df.families[,1] <- sub(pattern = "p.*",replacement ="p",x = df.families[,1])
df2 <- merge(df.families,df1,by=1)
df2 <- aggregate(df2[,-c(1:2)],by = df2["miR family"],FUN = sum,na.rm=T) #sum counts by family
gene_vector <- c("miR-150-5p","miR-25-3p/32-5p/92-3p/363-3p/367-3p","miR-142-3p.1",
"miR-17-5p/20-5p/93-5p/106-5p","miR-191-5p",
"miR-15-5p/16-5p/195-5p/322-5p/497-5p","miR-26-5p","miR-138-5p",
"miR-223-3p","miR-342-3p","miR-22-3p","miR-192-5p/215-5p",
"miR-125-5p/351-5p","miR-126-3p.1","miR-199-3p")
df3 <- subset(df2, df2$`miR family`%in% gene_vector) %>% arrange(factor(`miR family`, levels=gene_vector))
rownames(df3) <- df3[,1]
df3.scaled <- as.data.frame(t(scale(t(df3[-1])))) #Zscore normalization. Scaled only works on columns, so need to transform
setwd(fdir)
openxlsx::write.xlsx(df2.scaled,"Figure_1D_table.xlsx", rowNames=T)
#the df was exported and used to create a heatmap using Graphpad
CB2025_figure_3_RNAseq.R 0 → 100644
View file @ a78a5561
library(readxl)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(fgsea)
library(clusterProfiler)
library(enrichplot)
##### Directories #####
#us <- "/Users/Squadrito/"
us <-"C:/Users/bresesti.chiara/"
wdir<-paste0(us,"/Dropbox (HSR Global)/90-857433247_RNAseq_Squadrito/05-DGE-NoOut-Corr")
wdir_CB<-paste0(us, "/Dropbox (HSR Global)/90-857433247_RNAseq_Squadrito/Analysis CB_v2")
fdir <- paste0(us,"/Dropbox (HSR Global)/CancerGeneTherapy/Cancer Gene Therapy/MS/2024 Bresesti et al/Scripts/plots and tables used in figures")
#### Figure 3A&B ####
#Import df
setwd(wdir)
miR_ctrl <- read_excel("edgeR_results.xlsx", sheet = "miR342-control")
sponge_ctrl <- read_excel("edgeR_results.xlsx", sheet = "spongeBT-control")
#Add DEG color and label for volcano plot
miR_ctrl$DEG <- "NO"
miR_ctrl$DEG[miR_ctrl$logFC > 1 & miR_ctrl$PValue < 0.05] <- "UP"
miR_ctrl$DEG[miR_ctrl$logFC < (-1) & miR_ctrl$PValue < 0.05] <- "DOWN"
miR_ctrl$DEG_label <- NA
miR_ctrl$DEG_label[miR_ctrl$DEG != "NO"] <- miR_ctrl$...1[miR_ctrl$DEG != "NO"]
sponge_ctrl$DEG <- "NO"
sponge_ctrl$DEG[sponge_ctrl$logFC > 1 & sponge_ctrl$PValue < 0.05] <- "UP"
sponge_ctrl$DEG[sponge_ctrl$logFC < (-1) & sponge_ctrl$PValue < 0.05] <- "DOWN"
sponge_ctrl$DEG_label <- NA
sponge_ctrl$DEG_label[sponge_ctrl$DEG != "NO"] <- sponge_ctrl$...1[sponge_ctrl$DEG != "NO"]
#Volcano plot (left panel)
data <- miR_ctrl
ggplot(data=data, aes(x=logFC, y=-log10(PValue), col=DEG, label=DEG_label)) +
geom_point() +
theme_minimal() +
geom_text_repel() +
scale_color_manual(values=c("blue", "grey", "red")) +
geom_vline(xintercept=c(-1, 1), col="black") +
geom_hline(yintercept=-log10(0.05), col="black") +
scale_x_continuous(limits = c(-4.1,4.1), breaks = seq(-4,4,1)) +
scale_y_continuous(limits = c(0,20)) +
xlab(bquote(Log[2](FC))) +
ylab(bquote(-Log[10](PVal)))
data <- sponge_ctrl
ggplot(data=data, aes(x=logFC, y=-log10(PValue), col=DEG, label=DEG_label)) +
geom_point() +
theme_minimal() +
geom_text_repel() +
scale_color_manual(values=c("blue", "grey", "red")) +
geom_vline(xintercept=c(-1, 1), col="black") +
geom_hline(yintercept=-log10(0.05), col="black") +
scale_x_continuous(limits = c(-4.1,4.1), breaks = seq(-4,4,1)) +
scale_y_continuous(limits = c(0,20)) +
xlab(bquote(Log[2](FC))) +
ylab(bquote(-Log[10](PVal)))
#Import miR-342-3p target list (from TargetScan)
setwd(wdir_CB)
miR342_targets <- read_excel("TargetScan8.0__miR-342-3p.predicted_targets.xlsx")
miR342_targets <- filter(miR342_targets, miR342_targets$`Cumulative weighted context++ score`< (-0.3))
#ecdf plot (right panel)
data <- miR_ctrl
test <-filter(data, data$...1 %in% miR342_targets$'Target gene')
plot(ecdf(data$logFC), lwd = 2, do.points=F, verticals=T,
ylab = "Fraction of DEG",
xlab = bquote(Log[2](FC)),
xlim=c(-0.5,0.5),
ylim=c(0,1),
main="ECDF of gene expression: miR vs mut") +
lines(ecdf(test$logFC), col= "blue", do.points=F, lwd = 2, verticals=T) +
abline(v=0, col="black") +
abline(h=0.5, col="black") +
legend("bottomright", c("All genes","miR-342-3p targets"),
col = c("black","blue"), lwd=2, cex = 0.7)
ks.test(test$logFC, data$logFC, alternative = "g") #Kolmogorov-Smirnov test to calculate p-val of miR target distribution *not* greater than average data
data <- sponge_ctrl
test <-filter(data, data$...1 %in% miR342_targets$'Target gene')
plot(ecdf(data$logFC), lwd = 2, do.points=F, verticals=T,
ylab = "Fraction of DEG",
xlab = bquote(Log[2](FC)),
xlim=c(-0.3,0.3),
ylim=c(0,1),
main="ECDF of gene expression: sponge vs scr") +
lines(ecdf(test$logFC), col= "red", do.points=F, lwd = 2, verticals=T) +
abline(v=0, col="black") +
abline(h=0.5, col="black") +
legend("bottomright", c("All genes","miR-342-3p targets"),
col = c("black","red"), lwd=2, cex = 0.7)
ks.test(test$logFC, data$logFC, alternative = "l") #Kolmogorov-Smirnov test to calculate p-val of miR target distribution *not* less than average data
#### Figure 3C&D ####
#Upload df
setwd(wdir)
df1<- read_excel("EdgeR_results.xlsx",, sheet = "miR342-control")
setwd(wdir_CB)
miDB_sig5 <- readRDS("miDB_sig5.MLS.rds") #GO terms db
#Rename pathways of interest in miDB_sig5
names(miDB_sig5)[names(miDB_sig5) == "HALLMARK_OXIDATIVE_PHOSPHORYLATION"] <- "Oxydative phosphorylation"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_REGULATION_OF_CHOLESTEROL_METABOLIC_PROCESS"] <- "Regulation of cholesterol metabolic process"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_RESPONSE_TO_INTERLEUKIN_12"] <- "Response to IL12"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_REGULATION_OF_LIPID_BIOSYNTHETIC_PROCESS"] <- "Regulation of lipid biosynthetic process"
names(miDB_sig5)[names(miDB_sig5) == "GOMF_MHC_CLASS_I_PROTEIN_BINDING"] <- "MHC-I protein binding"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_TUMOR_NECROSIS_FACTOR_MEDIATED_SIGNALING_PATHWAY"] <- "TNFa mediated signalling pathway"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_NEGATIVE_REGULATION_OF_CELL_CYCLE_G2_M_PHASE_TRANSITION"] <- "Negative regulation of cell cycle progression"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_FATTY_ACYL_COA_BIOSYNTHETIC_PROCESS"] <- "Fatty acyl-CoA biosynthetic process"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_RESPONSE_TO_TUMOR_NECROSIS_FACTOR"] <- "Response to TNFa"
names(miDB_sig5)[names(miDB_sig5) == "HALLMARK_G2M_CHECKPOINT"] <- "G2M checkpoint"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_RESPONSE_TO_TYPE_I_INTERFERON"] <- "Response to type I interferon"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_INTERLEUKIN_1_MEDIATED_SIGNALING_PATHWAY"] <- "IL1 mediated signalling pathway"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_RESPONSE_TO_INTERLEUKIN_1"] <- "Response to IL1"
names(miDB_sig5)[names(miDB_sig5) == "LPS_RO"] <- "LPS response genes"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_CYTOKINE_MEDIATED_SIGNALING_PATHWAY"] <- "Cytokine mediated signalling pathway"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_MYELOID_CELL_DIFFERENTIATION"] <- "Myeloid cell differentiation"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_INTERLEUKIN_10_PRODUCTION"] <- "IL10 production"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_POSITIVE_REGULATION_OF_ENDOTHELIAL_CELL_PROLIFERATION"] <- "Regulation of endothelial cell proliferation"
names(miDB_sig5)[names(miDB_sig5) == "GOBP_TUMOR_NECROSIS_FACTOR_SUPERFAMILY_CYTOKINE_PRODUCTION"] <- "TNFa superfamily cytokine production"
names(miDB_sig5)[names(miDB_sig5) == "HALLMARK_ANGIOGENESIS"] <- "Angiogenesis"
names(miDB_sig5)[names(miDB_sig5) == "GOMF_PATTERN_RECOGNITION_RECEPTOR_ACTIVITY"] <- "Pattern recognition receptor activity"
names(miDB_sig5)[names(miDB_sig5) == "PGE2_RO"] <- "PGE2 response genes"
#Filter low expression, NA and order by FC
df1 <- filter(df1, df1$logCPM>5 & !is.na(df1$...1))
df1<-df1[order(df1$logFC),]
df2<-df1$logFC
names(df2)<-df1$...1
#Run GSEA with fgsea and filter by PVal
sig <- miDB_sig5
test1<-fgsea(sig, df2,minSize = 7,maxSize =500, nproc=1)
PvalResult<-filter(test1, test1$padj<= 0.05)
#DotPlot pathways of interest (Fig.3C)
pathways <- c("Oxydative phosphorylation",
"Regulation of cholesterol metabolic process",
"Response to IL12",
"Regulation of lipid biosynthetic process",
"MHC-I protein binding",
"TNFa mediated signalling pathway",
"Negative regulation of cell cycle progression",
"Fatty acyl-CoA biosynthetic process",
"Response to TNFa",
"G2M checkpoint",
"Response to type I interferon",
"IL1 mediated signalling pathway",
"Response to IL1",
"LPS response genes",
"Cytokine mediated signalling pathway",
"Myeloid cell differentiation",
"IL10 production",
"Regulation of endothelial cell proliferation",
"TNFa superfamily cytokine production",
"Angiogenesis",
"Pattern recognition receptor activity",
"PGE2 response genes")
genelist <- as.data.frame(PvalResult[PvalResult$pathway %in% pathways,])
ggplot(genelist, aes(x=NES, y=reorder(pathway,NES), size=size ,color=padj)) +
geom_point() +
scale_size_area(limits=c(10,450), max_size = 15) +
scale_colour_gradient(low="red",high="blue") +
labs(y='Pathway',x='NES')
#Run GSEA with ClusterProfiler
genelist <- data.frame(term = rep(names(miDB_sig5), sapply(miDB_sig5, length)),
gene = unlist(miDB_sig5))
df2 = sort(df2, decreasing = TRUE)
test2 <- GSEA(df2,TERM2GENE = genelist)
#Enrichment map (Fig.3D)
test3 <- filter(test2, ID %in% pathways)
test3 <- pairwise_termsim(test3)
emapplot(test3, min_edge = 0.01, color = "NES", layout="fr", repel =T)+
scale_fill_gradient2(name=bquote(NES),low="blue", high="red")
#N.B. Produces slightly different plot every time,
#but connection between pathways stays the same
\ No newline at end of file
CB2025_figure_5_scRNAseq.R 0 → 100644
View file @ a78a5561
library(Seurat)
library(openxlsx)
library(readxl)
library(ggplot2)
library(dplyr)
library(fgsea)
library(pheatmap)
library(scales)
library(RColorBrewer)
library(cetcolor)
#username <- "/Users/Squadrito/"
#username <- "/Users/bresesti.chiara/"
#wdir <- paste0(username, "Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/Analysis_MM/results2/CB1_CB3_CB4")
#wdir_CB <- paste0(username, "Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/CB_scripts/results2")
#plot_dir <- paste0(username,"/Dropbox (HSR Global)/CancerGeneTherapy/Cancer Gene Therapy/MS/2024 Bresesti et al/Scripts/plots and tables used in figures") # fdir previously
#dirCB <-paste0(username,"Dropbox (HSR Global)/90-935462466_scRNAseq_Bresesti/CB_scripts")
#sig <- readRDS(paste0(wdir, "/miDB_sig5.MLS.rds")) #GO terms db
wdir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/results2/CB1_CB3_CB4"
wdir_CB <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts/results2"
plot_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts_final/plots and tables used in figures"
dirCB <-"/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/CB_scripts"
dir.create(plot_dir, showWarnings=F, recursive=T)
# 16357 GO terms db. used for GSEA
sig <- readRDS("/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/reference/miDB_sig5.MLS.rds")
###############################################################################
set.seed(42)
obs1 <- readRDS(paste0(wdir, "/CB1_CB3_CB4_final.rds"))
# Annotate each cell with a sample label based on its original identity (library) and hash ID
obs1@meta.data$Sample <- "Undefined"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse1"] <- "Ctrl_01"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse2"] <- "Ctrl_02"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB1"&obs1@meta.data$hash.ID=="Mouse3"] <- "Ctrl_03"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse1"] <- "miR_01"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse2"] <- "miR_02"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB3"&obs1@meta.data$hash.ID=="Mouse3"] <- "miR_03"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse1"] <- "miR_04"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse2"] <- "miR_02"
obs1@meta.data$Sample[obs1@meta.data$orig.ident=="CB4"&obs1@meta.data$hash.ID=="Mouse3"] <- "miR_05"
obs1<-FindClusters(obs1, graph.name="RNA_snn_h.orig.ident", resolution=4)
# Manually annotate the clusters
obs1@meta.data$Annotation.final.CB<-"Undefined"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==0]<- "Monocytes"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==1]<- "TAMs_3"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==2]<- "TAMs_4"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==3]<- "LSECs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==4]<- "TAMs_1"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==5]<- "Cancer cells"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==6]<- "cDC2_MoDCs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==7]<- "TAMs_2"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==8]<- "T and NK cells"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==9]<- "Neutrophils"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==10]<- "pre DCs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==11]<- "KCs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==12]<- "TAMs_5"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==13]<- "B cells"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==14]<- "Hepatocytes"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==15]<- "Mig cDCs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==16]<- "Cholangiocytes"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==17]<- "cDC1"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==18]<- "Macrophages"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.0.6==19]<- "pDCs"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.1.2==8]<-"KC-LSEC doublet"
obs1@meta.data$Annotation.final.CB[obs1@meta.data$RNA_snn_h.orig.ident_res.4==51]<-"Basophils"
# Visualize final annotation in Umap (Fig.5A)
obs1 <- SetIdent(obs1,value="Annotation.final.CB")
pdf(file = paste(plot_dir, "CB_Annotation.final.CB_Fig.5A.pdf", sep = "/"), width=9, height=7)
DimPlot(obs1, reduction="umap.harmony.orig.ident", group.by="Annotation.final.CB", label=T, repel=T)
dev.off()
# Adding a simplified annotation
obs1@meta.data$Shortlables<-obs1@meta.data$Annotation.final.CB
DCs<-c("Mig cDCs","pDCs","pre DCs", "cDC1", "cDC2_MoDCs")
obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%DCs]<- "DCs"
granu<-c("Neutrophils","Basophils")
obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%granu]<- "Granulocytes"
TAMs<-c("TAMs_1", "TAMs_2", "TAMs_3", "TAMs_4", "TAMs_5")
obs1@meta.data$Shortlables[obs1@meta.data$Annotation.final.CB%in%TAMs]<- "TAMs"
# Dotplot with markers for population (Suppl.Fig.5A)
obs1 <- SetIdent(obs1,value="Shortlables")
obs1$Shortlables <- factor(x = obs1$Shortlables,
levels = c("Cancer cells","Cholangiocytes",
"Hepatocytes", "LSECs", "KC-LSEC doublet",
"KCs", "Macrophages","TAMs", "Monocytes",
"Granulocytes","DCs", "B cells", "T and NK cells"))
gene <- c("Vil1", "Cdx2","Cldn4","Cldn7","Epcam",
"Krt19", "Sox9", "Cftr", "Hnf1b",
"Alb", "Ttr", "Hnf4a", "Cyp3a11",
"Lyve1", "Stab2", "Fcgr2b", "Tek",
"Vsig4", "Clec4f", "Timd4",
"Cd5l", "Mertk", "Timd4", "Adgre1","Cd68","Marco",
"Itgam", "Cd14","Arg1", "Ccl2",
"Ccr2", "Itgam", "Cx3cr1",
"Ly6g", "S100a8", "S100a9","Mcpt8", "Cpa3",
"Cd209a","Itgax", "Flt3","H2-Aa","Xcr1", "Clec9a",
"Cd19", "Cd22", "Cd79a", "Ms4a1", "Pax5",
"Cd3e", "Cd4", "Cd8a", "Ncr1", "Gzmb")
pdf(file = paste(plot_dir, "CB_DotPlot_Shortlables_Fig.S5A.pdf", sep="/"), width=12, height=6)
DotPlot(obs1,features = unique(gene), cols = "RdBu") +
theme(axis.text.x = element_text(angle=45, vjust=1, hjust=1))
dev.off()
# Generate table with distribution clusters by Mouse
# Table saved in excel and used to generate a barplot with GraphPad Prism (Suppl.Fig.5B)
cells_labels <- table(obs1@meta.data$Annotation.final.CB, obs1@meta.data$Sample)
cells_labels_p <- round(prop.table(table(obs1@meta.data$Annotation.final.CB, obs1@meta.data$Sample), 2) * 100, 1)
write.csv(cells_labels, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_Fig.S5B.csv"), row.names=T)
write.csv(cells_labels_p, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_percentage_Fig.S5B.csv"), row.names=T)
#### Analysis of mOrange+ fraction (Fig.5B, Suppl.Fig.5C,D) ####
# The metadata indicating mOrange positivity are: positive_cells (generated with the mapping), mOrange_label (generated with CellRanger), mOrange_union (union of the previous two labels)
# Generate table with mOrange+ cells distribution by cluster
# Table saved in excel and used to generate a barplot with GraphPad Prism (Suppl.Fig.5C)
# Suppl.Fig.5D was subsequently generated dividing data from Suppl.Fig.5C and B
obs2 <- subset(obs1, subset = mOrange_union == "mOrange+")
cells_labels_orange <- table(obs2@meta.data$Annotation.final.CB, obs2@meta.data$Sample)
cells_labels_orange_p <- round(prop.table(table(obs2@meta.data$Annotation.final.CB, obs2@meta.data$Sample), 2) * 100, 1)
write.csv(cells_labels_orange, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_mOrange_fraction_Fig.S5C.csv"), row.names=T)
write.csv(cells_labels_orange_p, paste0(plot_dir, "/number_cells_sample_Annotation.final.CB_percentage_mOrange_fraction_Fig.S5C.csv"), row.names=T)
# DimPlot with overlayed mOrange density (Fig.5B)
data <- Embeddings(object = obs2[["umap.harmony.orig.ident"]])[(colnames(obs2)), c(1,2)]
data <- as.data.frame(data)
obs1 <- SetIdent(obs1, value="cells")
pdf(file = paste(plot_dir, "UMAP_mOrange_cells_Fig.5B.pdf", sep = "/"), width=9, height=7)
DimPlot(obs1, reduction = "umap.harmony.orig.ident", pt.size = 0.1, cols = "azure3") +
stat_density_2d(aes_string(x = "UMAPh_1", y = "UMAPh_2", fill = "after_stat(level)"), data=data,
linewidth = 0.4, geom = "density_2d_filled",
colour = "black", alpha = 0.2, n = 150, h = c(2,2))+
scale_fill_gradientn(colours = cet_pal(16, name = "fire"))
dev.off()
#### GSEA (Fig.5C-E) ####
# Subset populations of interest (to plot)
selected_pop <- c("DCs", "KCs", "TAMs", "Monocytes", "Cancer cells", "LSECs")
obs3 <- subset(obs1, subset = Shortlables %in% selected_pop)
# Dotplot of Slc7a11 expression by cluster and experimental group (Fig.5C, groups divided in Illustrator)
obs3 <- SetIdent(obs3, value="Shortlables")
obs3$Shortlables <- factor(x = obs3$Shortlables, levels = c("LSECs", "Cancer cells", "Monocytes", "TAMs", "KCs", "DCs"))
pdf(file = paste(plot_dir, "CB_DotPlot_Shortlables_Slc7a11_Fig.5C.pdf", sep = "/"), width=9, height=5)
DotPlot(obs3, features = "Slc7a11", dot.scale=20, cols = "RdBu", split.by = "PROP.Group") +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
coord_flip()
dev.off()
# Run DEG between groups in each pop
simpl_DEG <- createWorkbook()
deg_results <- list()
for(i in selected_pop){
obs3 <- SetIdent(obs3, value="Shortlables")
# We want to calculate DGE between miR342 vs miRmut in each cell population in Shortlables
markers <- FindMarkers(obs3, ident.1 = "miR342", group.by = "PROP.Group", subset.ident = i, min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3) # added min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3
markers$diff_pct <- abs(markers$pct.1-markers$pct.2)
name <-paste("DEG_simpl", i, sep = "_")
deg_results[[name]] <- markers
assign(name, markers)
addWorksheet(simpl_DEG, name)
writeData(simpl_DEG, name, markers, rowNames = T)
}
saveWorkbook(simpl_DEG, file = paste0(plot_dir, "/DEG_simpl_Shortlables.xlsx"), overwrite = T)
# Run GSEA for each population
GSEA_simpl <- data.frame(matrix(nrow = 0, ncol = 9))
for(i in selected_pop){
df1 <- deg_results[[paste0("DEG_simpl_", i)]]
df1 <- df1[order(df1$avg_log2FC),]
ranks<-df1$avg_log2FC
names(ranks)<-rownames(df1)
test<-fgsea(pathways = sig, stats = ranks, minSize=5, maxSize=1000, nproc=1, nPermSimple=10000)
test$cluster <- i
colnames(GSEA_simpl) <- colnames(test)
GSEA_simpl <- rbind(GSEA_simpl, test)
}
write.xlsx(GSEA_simpl, file = paste0(plot_dir, "/GSEA_simpl_Shortlables.xlsx"))
# Plot barplot for selected gene sets (Fig.5D, names and colors edited on illustrator)
gene_sets <- c("HALLMARK_TNFA_SIGNALING_VIA_NFKB",
"HALLMARK_INFLAMMATORY_RESPONSE",
"GOMF_CHEMOKINE_ACTIVITY",
"GOBP_RESPONSE_TO_CYTOKINE",
"GOBP_CYTOKINE_MEDIATED_SIGNALING_PATHWAY",
"HALLMARK_INTERFERON_GAMMA_RESPONSE",
"HALLMARK_INTERFERON_ALPHA_RESPONSE",
"GOBP_RESPONSE_TO_TYPE_I_INTERFERON",
"GOMF_ANTIOXIDANT_ACTIVITY",
"GOBP_CELLULAR_OXIDANT_DETOXIFICATION",
"GOMF_GLUTATHIONE_PEROXIDASE_ACTIVITY",
"GOMF_GLUTATHIONE_TRANSFERASE_ACTIVITY",
"GOBP_GLUTATHIONE_DERIVATIVE_METABOLIC_PROCESS",
"GOBP_GLUTATHIONE_METABOLIC_PROCESS",
"GOBP_ANTIGEN_PROCESSING_AND_PRESENTATION",
"GOMF_ANTIGEN_BINDING",
"GOMF_MHC_CLASS_II_PROTEIN_COMPLEX_BINDING",
"GOMF_MHC_PROTEIN_COMPLEX_BINDING",
"GOCC_MHC_CLASS_II_PROTEIN_COMPLEX",
"GOCC_MHC_PROTEIN_COMPLEX",
"GOBP_LEUKOCYTE_MIGRATION_INVOLVED_IN_INFLAMMATORY_RESPONSE",
"GOBP_CYTOKINESIS",
"GOBP_MICROTUBULE_CYTOSKELETON_ORGANIZATION",
"GOBP_MICROTUBULE_BASED_PROCESS",
"GOBP_POSITIVE_REGULATION_OF_CELL_CYCLE_PROCESS",
"GOBP_CELL_DIVISION",
"HALLMARK_G2M_CHECKPOINT",
"GOBP_MITOTIC_CELL_CYCLE",
"GOBP_CELL_CYCLE_PHASE_TRANSITION")
plot <- filter(GSEA_simpl, pathway %in% gene_sets)
plot$cluster <- factor(plot$cluster, levels = c("LSECs", "Cancer cells", "TAMs", "Monocytes","KCs","DCs"))
# Remove the 'HALLMARK_' prefix and other prefixes from the 'pathway' column
plot$pathway <- gsub(x=plot$pathway, pattern="HALLMARK_|GOBP_|GOMF_|GOCC_", replacement="")
plot$pathway <- gsub(x=plot$pathway, pattern="_", replacement=" ")
gene_sets <- gsub(x=gene_sets, pattern="HALLMARK_|GOBP_|GOMF_|GOCC_", replacement="")
gene_sets <- gsub(x=gene_sets, pattern="_", replacement=" ")
# Reorder the terms
plot$pathway <- factor(plot$pathway, levels = gene_sets)
write.xlsx(plot, file = paste0(plot_dir, "/GSEA_plotted.xlsx"))
# Plot the GSEA results (Fig5D)
p <- ggplot(data=plot, aes(x=NES, y=pathway, fill=log10(pval))) +
geom_bar(stat="identity", color="black")+
scale_fill_gradientn(
colors=c("red","white","blue"),
values=rescale(c(-8,log10(0.05),0)),
limits=c(-8,0),
na.value = "red") +
labs(y='Pathway', x='NES')+
theme(strip.text = element_text(face = "bold"),
axis.text.x = element_text(colour = "black"),
axis.text.y = element_text(colour = "black"))+
geom_vline(xintercept=0)+
facet_wrap(~cluster, ncol = length(unique(plot$cluster)))
ggsave(filename = paste0(plot_dir, "/GSEA_Fig5D.pdf"), plot=p, width=12, height=7)
# Create signature file with only signature cytokines for the more stringent GSEA in Fig.5E
sig.cyto <- sig[15061:16357]
# Define the populations we want to use for this GSEA
selected_pop2 <- c("B cells","cDC1","cDC2_MoDCs","KCs","Macrophages","Mig cDCs","Monocytes","Neutrophils","pDCs","T and NK cells", "TAMs_1","TAMs_2","TAMs_3","TAMs_4","TAMs_5")
# Run DEG between groups in each population
simpl_DEG2 <- createWorkbook()
deg_results <- list()
for(i in selected_pop2){
obs1 <- SetIdent(obs1, value="Annotation.final.CB")
# We want to calculate DGE between miR342 vs miRmut in each cell population in Annotation.final.CB
markers <- FindMarkers(obs1, ident.1="miR342", group.by="PROP.Group", subset.ident=i, min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3) # added min.pct=0.01, logfc.threshold=0, min.cells.feature=3, min.cells.group=3
markers$diff_pct <- abs(markers$pct.1-markers$pct.2)
name <-paste("DEG_simpl", i, sep = "_")
deg_results[[name]] <- markers
assign(name, markers)
addWorksheet(simpl_DEG2, name)
writeData(simpl_DEG2, name, markers, rowNames = T)
}
saveWorkbook(simpl_DEG2, file = paste0(plot_dir, "/DEG_Annotation.final.CB.xlsx"), overwrite = T)
# Repeat GSEA for cytokines signatures only
GSEA_cyto <- data.frame(matrix(nrow=0, ncol=9))
for(i in selected_pop2){
df1 <- deg_results[[paste0("DEG_simpl_", i)]]
df1<-df1[order(df1$avg_log2FC),]
ranks<-df1$avg_log2FC
names(ranks)<-rownames(df1)
test<-fgsea(pathways = sig.cyto, stats = ranks, minSize=2, maxSize=1000, nproc=1, nPermSimple=30000)
test$cluster <- i
colnames(GSEA_cyto) <- colnames(test)
GSEA_cyto <- rbind(GSEA_cyto,test)
}
write.xlsx(GSEA_cyto, file = paste0(plot_dir, "/GSEA_cyto.xlsx"))
# Subset the GSEA_cyto to match the cytokine signature calculated in the populations with the same populations that we have in our scRNAseq.
# This is necessary because each cytokine was tested and the effect in each population was saved as a different signature.
df <- GSEA_cyto
df_Bcells <- df[grepl("B_cell", df$pathway) & df$cluster=="B cells",]
df_cDC1 <- df[grepl("cDC1", df$pathway) & df$cluster=="cDC1",]
df_cDC2 <- df[grepl("cDC2", df$pathway) & df$cluster=="cDC2_MoDCs",]
df_macrophage <- df[grepl("Macrophage", df$pathway) & grepl("TAMs|KCs|Macrophages", df$cluster),]
df_MigDC <- df[grepl("MigDC", df$pathway) & df$cluster=="Mig cDCs",]
df_Monocyte <- df[grepl("Monocyte", df$pathway) & df$cluster=="Monocytes",]
df_Neutrophil <- df[grepl("Neutrophil", df$pathway) & df$cluster=="Neutrophils",]
df_pDC <- df[grepl("pDC", df$pathway) & df$cluster=="pDCs",]
df_Tcells <- df[grepl("T_cell_CD8", df$pathway) & df$cluster=="T and NK cells",]
# Combine again all the generated subsets
df <- rbind(df_Bcells, df_cDC1, df_cDC2, df_macrophage, df_MigDC, df_Monocyte, df_Neutrophil, df_pDC, df_Tcells)
# Remove the redundant data frames
rm(df_Bcells, df_cDC1, df_cDC2, df_macrophage, df_MigDC, df_Monocyte, df_Neutrophil, df_pDC, df_Tcells)
# Generate a column indicating the cytokine analysed.
df$cytokine <- gsub(".*_","",df$pathway)
df$padj[is.na(df$padj)] <- 1 # Convert the NA into padj=1
# Assing marks for the significance: *=padj<0.05; **=padj<0.01; ***=padj<0.001.
df$label <- ""
df <- df %>% mutate(label = cut(
padj,
breaks = c(0, 0.001, 0.01, 0.05, 1),
label = c("***", "**", "*", "")))
paletteLength <- 100
myColor <- colorRampPalette(rev(brewer.pal(5, "RdBu")))(paletteLength)
# Define breaks to set the zero to the center
df$NES[is.na(df$NES)] <- 0 # Convert the NA into NAS=0
MIN <- min(df$NES)
MAX <- max(df$NES)
myBreaks <- c(seq(MIN, 0, length.out=ceiling(paletteLength/2) + 1),
seq(MAX/paletteLength, MAX, length.out=floor(paletteLength/2)))
# Plot heatmap cytokine signatures (Fig.5E, removed dendrograms in illustrator)
mat.df <- with(df, tapply(NES, list(cluster, cytokine), FUN = mean))
mat.label <- with(df, tapply(label, list(cluster, cytokine), FUN = paste))
pheatmap(mat.df, scale = "none", display_numbers = mat.label,
clustering_method = "ward.D2",
cutree_rows = 2, cutree_cols = 3,
col=myColor, breaks = myBreaks, na_col = "grey",
filename = paste0(plot_dir, "/Heatmap_cytokine_sig_Fig.5E_new.pdf"), cellwidth=10, cellheight=10) # width=12, height=7
README.md 0 → 100644
View file @ a78a5561
# Bresesti_Cell_Reports_2025
**Reprogramming liver metastasis-associated macrophages towards an anti-tumoral phenotype through enforced miR-342 expression**
Bresesti C, Monti M, Beretta S, Merelli I, Squadrito M, _et al._; 2025
Corresponding Author: Mario Leonardo Squadrito. Email: squadrito.mario@hsr.it. [![Twitter URL](https://img.shields.io/twitter/url/https/twitter.com/wouter_decoster.svg?style=social&label=Follow%20%40Mariosqua)](https://x.com/Mariosqua)
GEO:
[GSE274043](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE274043) (single-cell RNA-seq)
[GSE274044](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE274044) (RNA-seq on iKCs)
[GSE274045](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE274045) (small RNA-seq on splenic and hepatic cell populations)
[GSE274046](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE274046) (bulk RNA-seq on splenic and hepatic cell populations)
## Directories and Files
- environment_singlecell5.yml: contains the conda virtual environment that can be used to install all the dependencies.
- scripts: folder with R scripts used for the analyses
- CB2025_figure_1_RNAseq.R
- CB2025_figure_3_RNAseq.R
- CB2025_figure_5_scRNAseq.R
- TCGA_analysis.R
- data: results of the analyses
- plots
- tables
- reference
- regev_lab_cell_cycle_mouse.rds: contains the genes involved in the cell cycle in mouse
- miDB_sig5.MLS.rds: reference files for GSEA analysis
- miDB_sig5.MLS_cytokines.rds: reference files for GSEA analysis
## scRNAseq analysis
The initial preprocessing of the data, including mapping against the _Mus musculus_ GRCm38 reference genome and gene counting, was done using the 10x Genomics Cell Ranger Software (v7.2.0) using default parameters. The resulting data were imported into R and analyzed with the Seurat package (v5.0.1).
TCGA_analysis.R 0 → 100644
View file @ a78a5561
library(data.table)
library(survival)
library(dplyr)
library(babelgene)
library(openxlsx)
library(readxl)
library(survminer)
library(gridExtra)
library(ggplot2)
##Load dataset
#username <- "C:/Users/notaro.marco/"
username <- "/Users/bresesti.chiara/"
wdir<- paste0(username, "/Dropbox (HSR Global)/CancerGeneTherapy/Cancer Gene Therapy/MS/2024 Bresesti et al Cell Reports/TCGA_analysis")
setwd(wdir)
metaD<-data.frame(fread("TCGA_phenotype_denseDataOnlyDownload.tsv.gz", sep='\t',colClasses=c("character"),data.table=FALSE))
Surv01<-data.frame(fread("Survival_SupplementalTable_S1_20171025_xena_sp", sep='\t',colClasses=c("character"),data.table=FALSE))
# FPKM01<-fread("tcga_RSEM_Hugo_norm_count")
FPKM01<-fread("pancanMiRs_EBadjOnProtocolPlatformWithoutRepsWithUnCorrectMiRs_08_04_16.xena.gz")
FPKM01<-as.data.frame(FPKM01)
# Annot<-data.frame(fread("probeMap_gencode.v23.annotation.gene.probemap", sep='\t',colClasses=c("character"),data.table=FALSE))
allGenes <- FPKM01[,1]
#######################################################################
#####PLOT HR Score of signature for multiple tumor types###############
#####################################################################
humanGenes<-list(c("hsa-miR-342-3p"))
####create xlsx file with genes, patients and tumor types
Pan.data3 <- metaD[,]
Pan.data4 <- merge(Surv01,Pan.data3,by=1)
Genes_selected <- FPKM01[FPKM01[,1]%in%humanGenes,]
rownames(Genes_selected)<-Genes_selected[,1]
Genes_selected<-Genes_selected[,-1]
Genes_selected <- data.frame(t(Genes_selected))
Genes_selected$sample<-rownames(Genes_selected)
Pan.data5 <- merge(Genes_selected,Pan.data4,by.y="sample")
colnames(Pan.data5)
#setwd(wdir)
#write.xlsx(Pan.data5, "miR342_patients_survival.xlsx")
#############Regression with defined k2
k2<-0.5 #define quartile for selecting low/high score
df2 <- NULL
for (i in unique(Pan.data5$cancer.type.abbreviation)){
Pan.data6 <- Pan.data5[Pan.data5$cancer.type.abbreviation%in%i,]
if(nrow(Pan.data6)>2){
P.OS <- Surv(as.numeric(Pan.data6$OS.time), as.numeric(Pan.data6$OS))
cohort.score<-rep(0,nrow(P.OS))
M1 <- quantile(Pan.data6$hsa.miR.342.3p,k2)
cohort.score[Pan.data6$hsa.miR.342.3p>=M1]<-1
a1<-coxph(P.OS ~ as.numeric(cohort.score))
df1 <- data.frame(Tumor = unique(Pan.data6$cancer.type.abbreviation),npatients=length(Pan.data6$X_PATIENT), HR = summary(a1)$coefficients[2], SE=summary(a1)$coefficients[3],
p = summary(a1)$coefficients[5] )
df2 <- rbind(df2,df1)
df2<-df2[(order(df2$p)),]
}else{}}
#write.xlsx(df2, "miR342_HR_by_tumortype.xlsx")
###Filters selected tumors
df3<-df2
df3<-df3[df3$p<0.1,]
df3<-df3[df3$npatients>100,]
df3$significance <- ifelse(df3$p < 0.05, "p<0.05", "NS")
a<-ggplot(df3, aes(x = reorder(Tumor, HR), y = HR, color = significance, label = significance)) +
geom_point(size = 4) +
geom_errorbar(aes(ymin = pmax(HR-SE, 0), ymax = HR + SE), width = 0.2) +
scale_color_manual(name = "Color", values = c("red", "blue")) +
labs(title = "HR by tumor type",
x = "Tumor Type",
y = "HR") +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 1), fill = "lightblue", alpha = 0.01) +
geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 1, ymax = Inf), fill = "pink", alpha = 0.01)+
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust= 1, color = "Black")) +
geom_hline(yintercept = 1, linetype = "dashed", color = "gray")
# pdf("HR_alltumors_342.pdf",width=9,height = 4)
a
# dev.off()
####Plot survival of chosen tumors
###Select type of tumor
df3<-df3[order(df3$HR),]
df3<-df3[df3$p<0.05,]
# pdf("Survival_***_342",width=10,height = 6)
par(mfrow = c(2,4))
survival<-list()
for (i in df3$Tumor){
tumor.type<-i
Pan.data6 <- Pan.data5[Pan.data5$cancer.type.abbreviation%in%i,]
P.OS <- Surv(as.numeric(Pan.data6$OS.time), as.numeric(Pan.data6$OS))
cohort.score<-rep(0,nrow(P.OS))
M1 <- quantile(Pan.data6$hsa.miR.342.3p,k2)
cohort.score[Pan.data6$hsa.miR.342.3p>=M1]<-1
##plotting survival
fit.score <- survfit(P.OS ~ cohort.score)
cox_model <- coxph(P.OS ~ as.numeric(cohort.score))
p_value <- summary(cox_model)$coefficients["as.numeric(cohort.score)", "Pr(>|z|)"]
formatted_p_value <- ifelse(p_value < 0.0001, "<0.0001", as.numeric(round(p_value,5)))
par(las = 0)
plot(fit.score, lty = c(1, 1), col = c("blue", "red"), xlab = "Time (d)", ylab = "Overall Survival", lwd = 2, bty = "n",
main = paste0("Survival in ", i))
# Add legend
legend(x = max(fit.score$time)/2, y = 1, legend = c("miRNA high", "miRNA low"), lty = c(1, 1), col = c("red", "blue"), lwd = 2,bty = "n")
text(x = max(fit.score$time)/1.8, y = 0.70, label = paste("p=", (formatted_p_value)), adj = c(0, 1))
text(x = 10, y = 0.05, label = paste("n =", length(cohort.score)), adj = c(0, 1))
survival[[i]] <- recordPlot()
}
# dev.off()
###Plot survival of metastastic/primary cancers
#k2<-0.5
Pan.data7 <- Pan.data5[Pan.data5$sample_type == "Metastatic",]
#Pan.data7 <- Pan.data5[Pan.data5$sample_type == "Primary Tumor",]
P.OS <- Surv(as.numeric(Pan.data7$OS.time), as.numeric(Pan.data7$OS))
cohort.score<-rep(0,nrow(P.OS))
M1 <- quantile(Pan.data7$hsa.miR.342.3p,k2)
cohort.score[Pan.data7$hsa.miR.342.3p>M1]<-1
fit.score <- survfit(P.OS ~ cohort.score)
cox_model <- coxph(P.OS ~ as.numeric(cohort.score))
p_value <- summary(cox_model)$coefficients["as.numeric(cohort.score)", "Pr(>|z|)"]
formatted_p_value <- ifelse(p_value < 0.0001, "<0.0001", as.numeric(round(p_value,5)))
par(las = 0)
pdf("Survival_metastatic_342.pdf",width=5,height =5)
plot(fit.score, lty = c(1, 1), col = c("blue", "red"), xlab = "Time (d)", ylab = "Overall Survival", lwd = 2, bty = "n",
main = "Metastatic tumors")
#Add legend
legend(x = max(fit.score$time)/2, y = 1, legend = c("miRNA high", "miRNA low"), lty = c(1, 1), col = c("red", "blue"), lwd = 2,bty = "n")
text(x = max(fit.score$time)/1.8, y = 0.70, label = paste("p=", (formatted_p_value)), adj = c(0, 1))
text(x = 10, y = 0.05, label = paste("n =", length(cohort.score)), adj = c(0, 1))
dev.off()
environment_singlecell5.yml 0 → 100644
View file @ a78a5561
This diff is collapsed.
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