library(Seurat)
library(dplyr)
library(tidyverse)
library(viridis)
library(ggalluvial)
library("ggsci")
library("ggplot2")
library("gridExtra")
library(ComplexHeatmap)
require(Matrix)
require(magrittr)
library(scales)
library(configr)
library(cowplot)
library(pheatmap)
library(RColorBrewer)
library(reshape2)
library(harmony)
library(SeuratData)
#################################
#Fig2
setwd('Fig2') #path
genes<-as.vector(read.table('Fig2B.genelist')$V1)
ccgg <- c('CLEC9A',
'WDFY4',
'IRF8',
'BATF3',
'IDO1',
'CLNK',
'CLEC10A',
'FCER1A',
'CD1C',
'CRIP1',
'CD207',
'S100B',
'CD1E',
'CD1A',
'HLA-DQB2',
'LAMP3',
'CD40',
'CD83',
'CD274',
'CCR7',
'MARCKSL1',
'IL4I1',
'TCF4',
'IRF4',
'PLAC8',
'SPIB',
'IL3RA')
DC<-readRDS('V1_without10_Celltype_allDC_2nd_withAnnotation.rds')
DC@meta.data$cell.type.sub3 <- "LAMP3+ DC"
DC@meta.data[DC@meta.data$cell.type.sub2 == "DC – CD1C+CLEC10A+", "cell.type.sub3"] <- "CD1c+ DC"
DC@meta.data[DC@meta.data$cell.type.sub2 == "DC – CD207+CD1E+", "cell.type.sub3"] <- "Langerin+ DC"
DC@meta.data[DC@meta.data$cell.type.sub2 == "DC – CLEC9A+WDFY4+", "cell.type.sub3"] <- "CD141+ DC"
DC@meta.data[DC@meta.data$cell.type.sub2 == "pDC – TCF4+IL3RA+", "cell.type.sub3"] <- "pDC"
DC@meta.data$cell.type.sub3 <- factor(DC@meta.data$cell.type.sub3,levels = c("LAMP3+ DC", "CD1c+ DC", "Langerin+ DC", "CD141+ DC", "pDC"))
celltypes<-c('CD141+ DC','CD1c+ DC','Langerin+ DC','LAMP3+ DC','pDC')
DC$cluster<-DC$cell.type.sub3
DC$cluster<-factor(DC$cluster,levels=c(celltypes))
Idents(DC)<-'cluster'
celllist<-c('CD141+ DC','CD1c+ DC','LAMP3+ DC','Langerin+ DC','pDC')
DC$cluster<-factor(DC$cluster,levels=celllist)
table(DC$cluster)
Idents(DC)<-'cluster'
DC$patient<-factor(DC$patient,levels=c('MIBC5','MIBC6','MIBC7','MIBC8','MIBC9','MIBC11','MIBC12','MIBC13','MIBC14','MIBC15','MIBC16'))
DC$status0<-DC$status
DC$status<- plyr::mapvalues(x =DC$status0,from = c('T','P'),to = c('malignant','non-malignant'))
DC$status<-factor(DC$status,levels=c('malignant','non-malignant'))
setdiff(genes,rownames(DC))
library(ComplexHeatmap)
library(circlize) #colorRamp2
colr<-c('blue','white','red')
col_fun = colorRamp2(c(-2,0,2), colr)
exprs <- t(data.frame(FetchData(object = DC, vars = genes)))
df1<- as.matrix(t(scale(t(exprs),center=T,scale=T)))# zscore ,默认对列进行scale
#use_raster=T
celltypes<-c('CD141+ DC','CD1c+ DC','Langerin+ DC','LAMP3+ DC','pDC')
cell1<-rownames(subset(DC@meta.data,cluster==celltypes[1]))
cell2<-rownames(subset(DC@meta.data,cluster==celltypes[2]))
cell3<-rownames(subset(DC@meta.data,cluster==celltypes[3]))
cell4<-rownames(subset(DC@meta.data,cluster==celltypes[4]))
cell5<-rownames(subset(DC@meta.data,cluster==celltypes[5]))
cells<-c(cell1,cell2,cell3,cell4,cell5)
meta<-DC@meta.data[cells,]
df<-df1[,cells]
column_ha<-HeatmapAnnotation(cluster=as.vector(meta$cluster),patient=as.vector(meta$patient),status=as.vector(meta$status))
# meta$status0<-meta$status
# meta$status<- plyr::mapvalues(x =meta$status0,from = c('T','P'),to = c('malignant','non-malignant'))
#celltype_colors
celltype_colors <- c("#FBD44A", "#8A9197","#709AE1", "#F07344", "#D2AF81")
names(celltype_colors)<-c('CD141+ DC','CD1c+ DC','Langerin+ DC','LAMP3+ DC','pDC')
#patient color
celllist<-c('CD141+ DC','CD1c+ DC','Langerin+ DC','LAMP3+ DC','pDC')
patientlist<-c('MIBC5','MIBC6','MIBC7','MIBC8','MIBC9','MIBC11','MIBC12','MIBC13','MIBC14','MIBC15','MIBC16')
library(RColorBrewer)
grid.col<-colorRampPalette(brewer.pal(8, "Set1"))(11)
names(grid.col)<-patientlist
status_colors<-c('#EE8084','#19BBC1')
names(status_colors)<-c('malignant','non-malignant')
ann_colors=list(cluster=celltype_colors,status=status_colors,patient=grid.col)
#column_ha
column_ha<-HeatmapAnnotation(cluster=factor(as.vector(meta$cluster),levels=celllist),patient=factor(as.vector(meta$patient),levels=patientlist),status=as.vector(meta$status),col = ann_colors,annotation_name_gp= gpar(fontsize = 8))
#Fig2B plot
pdf('Fig2B.DC_genes_ComplexHeatmap_new.pdf',w=10,h=8)
ComplexHeatmap::Heatmap(df,name = "Expression",cluster_columns = F,cluster_rows = F,use_raster=T,show_row_names=T,show_column_names=F,top_annotation = column_ha,col = col_fun)
dev.off()
#Fig2A plot
pdf("Fig2A.DC_umap_celltype.pdf",w=6.5,h=5)
p1 <- DimPlot(DC,raster=F, reduction = "umap_harmony",group.by = "cluster",pt.size = 1,cols=celltype_colors)
print(p1)
dev.off()
DC1<-subset(DC,cluster!='pDC')
DC1$cluster<-factor(DC1$cluster)
genes<-c('CD207','HLA-DRA','IDO1','CD86','CD80','CCR7','LAMP3','CLEC9A')
#Fig2D plot
p4<-VlnPlot(DC1,features=genes,pt.size=-1,group.by='cluster',slot="data",assay="RNA",stack = T, flip =T)+xlab('')+ggtitle('')+theme(legend.position='none',axis.text.x = element_text(color="black",size=12,angle=30))
pdf('Fig2D.DC_genes_Vlnplot_new.pdf',w=4,h=5)
print(p4)
dev.off()
#find marker genes
tmp0<-DC;prefix<-'DC_new'
Idents(tmp0)<-'cluster'
cluster.averages <- AverageExpression(object = tmp0, assays ='RNA',return.seurat = F)
write.csv(cluster.averages$RNA,file=paste(prefix,'all_averageExpression.csv',sep='_'),quote=F,row.names=T)
average.genes<-cluster.averages$RNA
Idents(tmp0)<-'cluster'
all.markers<-FindAllMarkers(tmp0, only.pos = T, min.pct =0.1, logfc.threshold = 0.25,test.use='wilcox')
markers<-subset(all.markers,p_val_adj<0.01)
write.csv(markers,file=paste0(prefix,'_all.markers_FDR0.01.csv'),quote=F,row.names=T)
#Fig2H monocle2 plot
require(monocle)
require(Seurat)
require(dplyr)
require(Matrix)
require(magrittr)
library(scales)
library(ggplot2)
library(configr)
library(cowplot)
HSMM_myo<-readRDS('DC3_HSMM_myo.rds')
#legend.title=element_blank(),
p0<-theme(legend.position='right',
axis.text.x = element_text(color="black",size=10),
axis.text.y = element_text(color="black",size=10),
axis.title.x = element_text(face="plain", color="black",size=15),
axis.title.y = element_text(face="plain", color="black",size=15))
prefix<-'DC3_new'
cluster_cols<-c("#FEE500","#8A9FD1","#C06CAB")
p3<-plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 1,show_branch_points = F,show_backbone = T)+p0+theme(legend.position='right')+ scale_colour_continuous(low='#1D068D',high ='#F4EA26', name = 'Pseudotime',guide=guide_colorbar(reverse=F))
graph <- paste(prefix,'plot_cell_trajectory_Pseudotime_new.pdf',sep='_')
pdf(graph, w=6, h=4)
print(p3)
dev.off()
celltype_colors <- c("#FBD44A", "#8A9197","#709AE1", "#F07344", "#D2AF81")
names(celltype_colors)<-c('CD141+ DC','CD1c+ DC','Langerin+ DC','LAMP3+ DC','pDC')
graph <- paste(prefix,'plot_cell_trajectory_celltype_new.pdf',sep='_')
pdf(graph, w=6.5, h=4)
p1<-plot_cell_trajectory(HSMM_myo, color_by = "cluster",cell_size=1,show_branch_points = F,show_backbone = T)+scale_color_manual(values=celltype_colors)+theme(text = element_text(size=15))+p0
print(p1)
dev.off()
#################################
#FigS2A
pdf("FigS2A.DC_umap_status.pdf",w=6.5,h=5)
p1 <- DimPlot(DC,raster=F, reduction = "umap_harmony",group.by = "status",pt.size = 1)
print(p1)
dev.off()
#################################
#Fig3A
setwd('Fig3')
#T umap
# top marker genes
Trds<-readRDS('V1_without10_Celltype_T_clean_5thAnnotation.rds')
Trds$cluster<-Trds$cell.type.sub
celltypes1<-c("naïve","naïve-RPL","CD4-Th1-like","CD4-Th17-1","CD4-Th17-2","CD4-Treg(naïve)","CD4-Treg(Th17)","CD8-TRM-1","CD8-TEM","CD8_TEMRA","CD8-CT-1","CD8-CT-2","CD8-IFNsignaling","CD8-IEL_like","CD8-UK","gdT-antitumor")
Trds$cluster<-factor(Trds$cluster,levels=celltypes1)
library(RColorBrewer)
cluster_cols<-colorRampPalette(brewer.pal(8, "Set1"))(16)
table(Trds$cluster)
Idents(Trds)<-'cluster'
pdf("Fig3A.T_umap_celltype_new1.pdf",w=6.5,h=5)
p1 <- DimPlot(Trds,raster=F, reduction = "umap",group.by = "cluster",pt.size = 0,cols=cluster_cols)
print(p1)
dev.off()
Trds$status0<-Trds$status
Trds$status<- plyr::mapvalues(x =Trds$status0,from = c('T','P'),to = c('malignant','non-malignant'))
Trds$status<-factor(Trds$status,levels=c('malignant','non-malignant'))
pdf("FigS3A.T_umap_status_new.pdf",w=6.5,h=5)
p1 <- DimPlot(Trds,raster=F, reduction = "umap_harmony",group.by = "status",pt.size = 0)
print(p1)
dev.off()
library(RColorBrewer)
cluster_cols<-colorRampPalette(brewer.pal(8, "Set1"))(16)
genes<-as.vector(read.table('FigS3B.gene.list')$V1)
ccgg2 <- ('SELL',
'CCR7',
'IL2RA',
'TNFRSF4',
'FOXP3',
'CD69',
'GZMA',
'GZMB',
'GZMK',
'PRF1',
'IL17A',
'CXCL13',
'HSPA1A',
'TIGIT',
'MKI67',
'STMN1')
p4<-VlnPlot(Trds,features=genes,pt.size=-1,group.by='cluster',slot="data",assay="RNA",stack = T, flip =T)+xlab('')+ggtitle('')+theme(legend.position='none',axis.text.x = element_text(color="black",size=8,angle=30))
pdf('FigS3B.T_genes_Vlnplot_new.pdf',w=6,h=7)
print(p4)
dev.off()
#############
#plot_sank
library('ggalluvial')
plot_sank <- function(dada_per, condition, groups, count, colors){
p0<-theme(panel.grid=element_blank(), legend.background = element_rect(colour = NA),
legend.title = element_blank(),legend.text = element_text(face="plain", color="black",size=8),
axis.text.x = element_text(color="black",size=12),
axis.text.y = element_text(color="black",size=10),
axis.title.x = element_text(face="plain", color="black",size=15),
axis.title.y = element_text(face="plain", color="black",size=15))
p <- ggplot(dada_per,
aes_string(x = condition, stratum = groups, alluvium = groups, y=count,
fill = groups, label = groups)) +
scale_x_discrete(expand = c(0, 0)) +
geom_flow(width = 1/8) +
geom_stratum(alpha = .9,width = 1/10) +
#geom_text(stat = "stratum", size = 3, color="black") +
scale_fill_manual(values = colors) +
xlab("") + ylab("") +
theme_bw() +
theme(panel.grid =element_blank()) +
theme(panel.border = element_blank()) +
theme(legend.title = element_blank()) +
theme(axis.line = element_blank(),axis.ticks = element_blank()) +p0 + ggtitle("")
return(p)
}
#Fig3B
library(RColorBrewer)
cluster_cols<-colorRampPalette(brewer.pal(8, "Set1"))(16)
meta<-Trds0@meta.data
per <- round(prop.table(table(meta[,c("status", "cluster")]), margin=1),3)
per <- data.frame(per)
plot_sank(per, "status", "cluster","Freq", cluster_cols)
ggsave('Fig3B.T_celltype_ratio_group_plot_sank_new.pdf',w=5,h=6)
#Fig3C
#find marker genes
tmp0<-Trds;prefix<-'T_new'
Idents(tmp0)<-'cluster'
cluster.averages <- AverageExpression(object = tmp0, assays ='RNA',return.seurat = F)
write.csv(cluster.averages$RNA,file=paste(prefix,'all_averageExpression.csv',sep='_'),quote=F,row.names=T)
average.genes<-cluster.averages$RNA
Idents(tmp0)<-'cluster'
all.markers<-FindAllMarkers(tmp0, only.pos = T, min.pct =0.1, logfc.threshold = 0.25,test.use='wilcox')
markers<-subset(all.markers,p_val_adj<0.01)
write.csv(markers,file=paste0(prefix,'_all.markers_FDR0.01.csv'),quote=F,row.names=T)
library(pheatmap)
topn<-5
top_gene <- markers %>% group_by(cluster) %>% top_n(topn,avg_log2FC)
filename<-paste0(prefix,'_top5_markers_genes_pheatmap.pdf')
genes<-unique(as.vector(top_gene$gene))
data<-average.genes
pheatmap(log2(data[genes,]+1),filename=filename,width=8,height=12,cluster_rows=F,cluster_cols=F,display_numbers = F,number_format = "%.0f",fontsize = 8,fontsize_col = 10,border_color=NA,angle_col = "45",scale='row')
#################################
#Fig6A
#MIBC5、MIBC6、MIBC8、MIBC13
DC$treat <- ifelse(DC$patient %in% c('MIBC5','MIBC6','MIBC8','MIBC13'),'Chemtherpy','No Chemtherpy')
pdf("Fig6A.DC_umap_treat_new.pdf",w=6.5,h=5)
p1 <- DimPlot(DC,raster=F, reduction = "umap_harmony",group.by = "treat",pt.size = 1)+ggtitle('DC')
print(p1)
dev.off()
#Fig6B
#MIBC5、MIBC6、MIBC8、MIBC13
Trds$treat <- ifelse(Trds$patient %in% c('MIBC5','MIBC6','MIBC8','MIBC13'),'Chemtherpy','No Chemtherpy')
pdf("Fig6B.T_umap_treat_new.pdf",w=6.5,h=5)
p1 <- DimPlot(Trds,raster=F, reduction = "umap_harmony",group.by = "treat",pt.size = 0)+ggtitle('T cells')
print(p1)
dev.off()
