Add files via upload

Figure1.R

@@ -0,0 +1,91 @@
+options(stringsAsFactors = F)
+library(cowplot)
+library(ggpubr)
+library(dplyr)
+library(readr)
+library(tidyr)
+library(ggforce)
+library(pals)
+library(pheatmap)
+library(scales)
+library(ggthemes)
+library(Seurat)
+
+## Figure 1
+
+## Fig1b: inflammation-associated genes
+inflammationScore <- function(x, cols=NULL){
+  inflgenes <- unlist(strsplit("IFNG, IFNGR1, IFNGR2, IL10, IL12A, IL12B, IL12RB1, IL12RB2, IL13, IL17A, IL17F, IL18, IL18R1, IL18RAP, IL1A, IL1B, IL2, IL21, IL21R, IL22, IL23A, IL23R, IL2RG, IL4, IL4R, IL5, IL6, JUN, NFKB1, RELA, RORA, RORC, S100A8, S100A9, STAT1, STAT3, STAT4, STAT6, TGFB1, TGFB2, TGFB3, TNF", ", *"))
+  Idents(x) <- x$sample
+  groupExp <- AverageExpression(x, assays = 'SCT', slot = 'data')$SCT
+  inflgenes <- inflgenes[inflgenes %in% rownames(groupExp)]
+  gMeanExp <- groupExp[inflgenes, ]
+  gMeanExp <- t(apply(gMeanExp, 1, function(x){(x - min(x)) / max(x)}))
+  gMeanExp <- reshape2::melt(t(gMeanExp)) %>% 
+    dplyr::rename(Sample="Var1", Gene="Var2")
+  p <- ggboxplot(gMeanExp, x = "Sample", y = "value", color = "Sample", 
+                 orientation = "horizontal", add = "jitter", 
+                 ylab = "Scaled mean", palette = cols) 
+  p
+}
+sampleCols <- setNames(c("#00AFBB", "#E7B800", "#E7B800", "#E7B800", "#FC4E07", "#FC4E07", "#FC4E07", "#FC4E07"),
+                       c("NT-P2", "PT-P1", "PT-P2", "PT-P3", "HM-P1", "HM-P2", "HM-P3", "HM-P4"))
+fig1b <- inflammationScore(seuratObj, sampleCols)
+fig1b
+
+## Fig1c: 29 clusters
+clusterCols <- c("#843C39", "#8C6D31", "#E6550D", "#3182BD", "#54990F", "#BD9E39", "#E7BA52", "#31A354", "#E41A1C", 
+                          "#6BAED6", "#9ECAE1", "#AD494A", "#E7CB94", "#74C476", "#A1D99B", "#C7E9C0", "#99600F", 
+                          "#E7298A", "#C3BC3F", "#D6616B", "#FF7F00", "#1B9E77", "#FDAE6B", "#B3823E", "#66A61E", 
+                          "#F1788D", "#C6DBEF", "#E6550D", "#E7969C")
+setNames(clusterCols) <- sprintf("%02d",0:28)
+fig1c <- DimPlot(seuratObj, group.by='cluster', cols=clusterCols, pt.size=1, raster=F) 
+
+
+## Fig1d 
+Idents(seuratObj) <- seuratObj$cluster
+avg <- AverageExpression(seuratObj, features=mmarkers, assays='SCT')$SCT 
+sdat <- t(apply(avg, 1, scale))
+colnames(sdat) <- colnames(avg)
+phtm <- pheatmap::pheatmap(sdat, cluster_rows = F, silent = T)
+cord <- phtm$tree_col$labels[phtm$tree_col$order]
+cord <- data.frame(cluster=factor(PHmeta, levels = unique(PTmeta[PTorder])), 
+                   id=factor(names(PHmeta),levels = cord)) %>% 
+  arrange(cluster, id) %>% pull(id) %>% as.character()
+geneOrder <- data.frame(cluster=factor(colnames(markerExp)[apply(markerExp, 1, which.max)], levels = cord))
+geneOrder$gene <- unlist(lapply(levels(geneOrder$cluster), function(i){
+  order(markerExp[which(as.character(geneOrder$cluster)==i),i], decreasing = T)
+}))
+markerExp <- markerExp[order(geneOrder$cluster, geneOrder$gene), ]
+
+mxid <- factor(colnames(sdat)[apply(sdat, 1, which.max)], levels=colnames(sdat))
+gord <- rownames(sdat)[order(mxid)]
+gord <- intersect(rownames(markerExp), gord)
+Idents(seuratObj) <- factor(seuratObj$cluster, levels = cord)
+fig1d <- DotPlot(seuratObj, features = gord) + rotate() +
+  scale_color_gradientn(colours=rev(brewer.rdylbu(20)), guide = "colourbar") +
+  theme_minimal_grid() + theme(axis.text.x=element_text(angle=90, hjust=1))
+fig1d
+
+## Fig1e 
+groupCols <- c(NT="#00AFBB", PT="#E7B800", HM="#FC4E07" )
+cellstat <- table(seuratObj$sample, seuratObj$majorcluster)
+cellstat <- as.data.frame(cellstat * 100 / rowSums(cellstat))
+colnames(cellstat) <- c("sample", "cluster", "percentage")
+cellstat$group <- factor(sampleGroups[as.character(cellstat$sample)], levels=names(groupCols))
+fig1e <- ggbarplot(cellstat, x = "cluster", y = "percentage", fill = "group",
+                   color = 'black', add = "mean_se", palette = groupCols,
+                   position = position_dodge(0.8)) + rotate()
+
+## Fig1g 
+majorColors <- setNames(c("#843C39", "#8C6D31", "#E6550D", "#3182BD", "#54990F", "#E41A1C", "#E7298A", "#C3BC3F", "#FF7F00", "#1B9E77", "#66A61E", "#F1788D"), 
+                        c("Ductal cell", "T cell","Fibroblasts", "NKT", "MKI67+ ductal cell",
+                          "Mast cell","Endocrine cell", "B cell","Plasma cell","Endothelial cell"))
+fig1g1 <- DimPlot(seuratObj, group.by='majorcluster', cols=majorColors, pt.size=1, raster=F) 
+fig1g1
+
+fig1g2 <- FeaturePlot(seuratObj, features='DEgenes', reduction='umap', order=TRUE, 
+                     cols=brewer.ylorrd(max(seuratObj$DEgenes)), 
+                     min.cutoff='q1', pt.size=1, combine=T, raster = F)
+fig1g2
+

Figure2.R

@@ -0,0 +1,84 @@
+## Figure 2
+
+## "Ductal cell", "MKI67+ ductal cell"
+epiObj <- subset(x = seuratObj, subset = majorcluster %in% c("Ductal cell","MKI67+ ductal cell"))
+epiObj <- RunUMAP(epiObj, verbose=TRUE) %>% 
+  FindNeighbors() %>% 
+  FindClusters(resolution = 0.25)
+DefaultAssay(epiObj) <- "SCT"
+
+
+## Fig2a 
+ecols <- setNames(c("#1B9E77", "#D95F02", "#7570B3", "#E7298A", "#66A61E", "#E6AB02"), c("0", 1:5))
+mkcols <- c("#222222", "#F3C300", "#875692", "#F38400")
+
+epiObj$flag <- ifelse(epiObj@assays$SCT@data['EPCAM',] > 0 & 
+                         epiObj@assays$SCT@data['MKI67',] > 0, 'EPCAM+MKI67+', 
+                       ifelse(epiObj@assays$SCT@data['EPCAM',] > 0, 'EPCAM+',
+                              ifelse(epiObj@assays$SCT@data['MKI67',] > 0, 'MKI67+','Other')))
+epiObj$flag <- factor(epiObj$flag, levels = c("Other",'MKI67+','EPCAM+MKI67+','EPCAM+'))
+
+cell_eb <- data.frame(epiObj@reductions[["umap"]]@cell.embeddings[,1:2], cluster=as.character(epiObj$cluster),
+                      marker=epiObj$flag, group=as.character(epiObj$group), stringsAsFactors = F)
+colnames(cell_eb) <- c('x','y','cluster','marker','group')
+pp1 <- ggplot(data=cell_eb, aes(x,y)) + geom_point(aes(colour = factor(marker)), size=0.25) +
+  scale_color_manual(values = mkcols) + theme_pubr() +
+  NoAxes() + NoLegend()
+pp2 <- ggplot(data=cell_eb, aes(x,y)) + geom_point(aes(colour = factor(cluster)), size=0.25) +
+  scale_color_manual(values = ecols) + theme_pubr() +
+  NoAxes() + NoLegend()
+
+cell_eb <- data.frame(epiObj@reductions[["umap"]]@cell.embeddings[,1:2], cluster=as.character(epiObj$cluster),
+                      group=as.character(epiObj$group), stringsAsFactors = F)
+colnames(cell_eb) <- c('x','y','cluster', 'group')
+pp3 <- ggplot(data=cell_eb, aes(x,y)) + geom_point(aes(colour = factor(group)), alpha=0.5, size=0.25) +
+  scale_color_manual(values = alpha(groupCols, alpha = 0.3)) + theme_pubr() +
+  NoAxes() + NoLegend()
+df <- epiObj@meta.data %>% dplyr::count(cluster, group) %>% dplyr::rename(cluster = cluster)
+df$n <- df$n/as.numeric(table(epiObj@meta.data$group)[df$group]) # normalize 组的细胞数量
+coords <- cell_eb %>% group_by(cluster) %>% dplyr::summarise(x=median(x), y=median(y), .groups = "keep") %>% as.data.frame
+coords <- coords %>% mutate(x=ifelse(cluster=='6',x+1,x), y=ifelse(cluster=='6',y+1,y))
+piedf <- reshape2::acast(df, cluster ~ group, value.var = "n")
+piedf[is.na(piedf)] <- 0
+piedf <- cbind(piedf, coords[match(rownames(piedf), coords$cluster), -1])
+cellnumbers <- dplyr::count(cell_eb, cluster)
+sizee <- cellnumbers[match(rownames(piedf), cellnumbers[,1]), 2]
+piedf$radius <- scales::rescale(sizee, c(0.4,0.9))
+pp3 <- pp3 + ggnewscale::new_scale_fill() + 
+  geom_scatterpie(aes_(x=~x,y=~y,r=~radius), data=piedf,
+                  cols=colnames(piedf)[1:(ncol(piedf)-3)],color=NA) +
+  coord_equal() + scale_fill_manual(values = groupCols)
+
+
+cell_eb <- data.frame(epiObj@reductions[["umap"]]@cell.embeddings[,1:2], cluster=as.character(epiObj$cluster),
+                      group=as.character(epiObj$CNV), stringsAsFactors = F)
+colnames(cell_eb) <- c('x','y','cluster', 'group')
+pp4 <- ggplot(data=cell_eb, aes(x,y)) + geom_point(aes(colour = factor(group)), alpha=0.5, size=0.25) +
+  scale_color_manual(values = alpha(cnvCols, alpha = 0.3)) + theme_pubr() +
+  NoAxes() + NoLegend()
+df <- epiObj@meta.data %>% dplyr::count(cluster, CNV) %>% dplyr::rename(cluster = cluster)
+df$n <- df$n/as.numeric(table(epiObj@meta.data$CNV)[df$CNV]) # normalize 组的细胞数量
+coords <- cell_eb %>% group_by(cluster) %>% dplyr::summarise(x=median(x), y=median(y), .groups = "keep") %>% as.data.frame
+coords <- coords %>% mutate(x=ifelse(cluster=='6',x+1,x), y=ifelse(cluster=='6',y+1,y))
+piedf <- reshape2::acast(df, cluster ~ CNV, value.var = "n")
+piedf[is.na(piedf)] <- 0
+piedf <- cbind(piedf, coords[match(rownames(piedf), coords$cluster), -1])
+cellnumbers <- dplyr::count(cell_eb, cluster)
+sizee <- cellnumbers[match(rownames(piedf), cellnumbers[,1]), 2]
+piedf$radius <- scales::rescale(sizee, c(0.4,0.9))
+pp4 <- pp4 + ggnewscale::new_scale_fill() + 
+  geom_scatterpie(aes_(x=~x,y=~y,r=~radius), data=piedf,
+                  cols=colnames(piedf)[1:(ncol(piedf)-3)],color=NA) +
+  coord_equal() + scale_fill_manual(values = cnvCols)
+
+
+print(ggarrange(as_ggplot(get_legend(p1)), as_ggplot(get_legend(p2)), as_ggplot(get_legend(p3)), as_ggplot(get_legend(p4)), ncol=2, nrow=2))
+print(ggarrange(pp1+NoLegend()+NoAxes()+ggtitle(NULL), 
+                pp2+NoLegend()+NoAxes()+ggtitle(NULL),
+                ncol=2, nrow=2))
+print(ggarrange(pp3+NoLegend()+NoAxes()+ggtitle(NULL), 
+                pp4+NoLegend()+NoAxes()+ggtitle(NULL), 
+                ncol=2, nrow=2))
+
+## Fig2b: epi_split_p_heatmap
+

Figure3.R

@@ -0,0 +1,84 @@
+## Figure 3
+
+## Fig3a 
+Idents(epiObj) <- factor(epiObj$cluster, levels = c("5",4,0,1,2,3))
+fig3a <- VlnPlot(epiObj, features=c("MKI67","EPCAM","CEACAM6","CEACAM5","KLK7"), 
+             stack = T, fill.by='ident', cols = ecols, pt.size=0) +
+  NoLegend() 
+print(fig3a)
+
+## Fig3b
+# see Velocity.ipynb as an example 
+
+## Fig3c 
+p1 <- monocle::plot_cell_trajectory(cds, cell_size=0.5, color_by="Pseudotime", show_branch_points=F) +
+  gradient_color(kovesi.rainbow_bgyr_35_85_c73(100))
+p2 <- monocle::plot_cell_trajectory(cds, cell_size=0.1, color_by="cluster", show_branch_points=F) + 
+  scale_color_manual(values=cols)
+p3 <- monocle::plot_cell_trajectory(cds, cell_size=0.1, color_by="group", show_branch_points=F) + 
+  scale_color_manual(values=groupCols)
+p4 <- monocle::plot_cell_trajectory(cds, cell_size=0.5, 
+                                    color_by="cluster", show_branch_points=F) + 
+  scale_color_manual(values=cols) + facet_wrap(~group)
+
+fig3c <- ggarrange(p1 + NoAxes() + NoLegend(), 
+                   p2 + NoAxes() + NoLegend(), 
+                   p3 + NoAxes() + NoLegend(),
+                   monocle::plot_cell_trajectory(cds, cell_size=0.1, color_by="State", show_branch_points=F) + 
+                     scale_color_manual(values=stateColor) + NoAxes() + NoLegend(), 
+                   monocle::plot_cell_trajectory(cds, cell_size=0.1, color_by="CNV", show_branch_points=F) + 
+                     scale_color_manual(values=cnvCols) + NoAxes() + NoLegend(), 
+                   plot_expression_trajectory('KLK7') + 
+                     gradient_color(c('grey',brewer.reds(10))) + NoAxes() + NoLegend(), 
+                   nrow=2, ncol=3)
+fig3c
+
+
+## Fig3e 
+cvnCols <- setNames(c("#b9ca5d", "#00a2b3", "#f3a546", "#cf3e53"), 
+                    c("Normal", "Low", "Medium", "High"))
+
+pdf("Fig3e1.pdf", height=4, width=8.27)
+op <- par(mar=rep(0.5,4), mfrow=c(4,1))
+x <- table(pltd$index, pltd$cluster)
+x <- x / rowSums(x)
+rownames(x) <- NULL
+barplot(t(x), col=ecols, border = NA, axes = F, main="Cluster")
+
+x <- table(pltd$index, pltd$group)
+x <- x / rowSums(x)
+rownames(x) <- NULL
+barplot(t(x), col=groupCols, border = NA, axes = F, main="Group")
+
+pltd$CNV <- epiObj@meta.data[rownames(pData(cds)),'CNV']
+x <- table(pltd$index, pltd$CNV)
+x <- x / rowSums(x)
+rownames(x) <- NULL
+barplot(t(x), col=cnvCols, border = NA, axes = F, main="CNV")
+
+x <- table(pltd$index, pltd$State)
+x <- x / rowSums(x)
+rownames(x) <- NULL
+barplot(t(x), col=stateColor, border = NA, axes = F, main="State")
+par(op)
+dev.off()
+
+epiObj@meta.data[,'bin'] <- pData(cds)[Cells(epiObj),'index']
+Idents(epiObj) <- 'bin'
+bindata <- AverageExpression(epiObj, assays = 'SCT', features=genex)$SCT %>% as.matrix()
+avg <- t(apply(bindata, 1, function(x){x=log1p(x);x/max(x)}))
+colnames(avg) <- NULL
+pdf("Fig3e2.pdf", height=8.27, width=8.27)
+pheatmap(avg, color = mapal, cluster_cols = F, cluster_rows = F, border_color = NA)
+dev.off()
+
+
+## Fig3f 
+fig3f <- monocle::plot_cell_trajectory(cds, cell_size=0.5, color_by="group") +
+  scale_color_manual(values=groupCols)
+fig3f
+
+
+## fig3g
+# see Velocity.ipynb
+

Figure4.R

@@ -0,0 +1,23 @@
+## Figure 4
+
+## Fig4a 
+fig4a <- DimPlot(cafObj, pt.size = 1, cols = fcols, label = F) + 
+  NoAxes() + NoLegend() + ggtitle(NULL)
+Fig4a
+
+
+## Fig4b 
+# similar to Fig3a and Fig1d
+
+## Fig4c 
+# similar to Fig2a
+
+## Fig4d 
+# similar to Fig3c
+
+## Fig4e 
+# similar to Fig1d
+
+## Fig4g-i 
+# similar to Fig3 
+

Figure5.R

@@ -0,0 +1,2 @@
+## Figure 5
+# refer to Figures 1-3

Figure6.R

@@ -0,0 +1,2 @@
+## Figure 6
+# refer to Figures 1-4