介绍

文章提供画图代码和数据，本文记录

数据和代码

数据可从以下链接下载（画图所需要的所有数据）：

• 百度云盘链接: https://pan.baidu.com/s/1peU1f8_TG2kUKXftkpYqug

• 提取码: 7pjy

图1

#### Figure 1: Census of cell types of the mouse uterine tube ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)#### Distal and Proximal Datasets ####Distal <- readRDS(file = "../dataset/Distal_Filtered_Cells.rds" , refhook =  NULL)Proximal <- readRDS( file = "../dataset/Proximal_Filtered_Cells.rds" , refhook =  NULL)#### Figure 1b: Cells of the Distal Uterine Tube ####Distal_Named <- RenameIdents(Distal, '0' = "Fibroblast 1", '1' = "Fibroblast 2", '2' = "Secretory Epithelial",'3' = "Smooth Muscle", '4' = "Ciliated Epithelial 1", '5' = "Fibroblast 3", '6' = "Stem-like Epithelial 1",'7' = "Stem-like Epithelial 2",'8' = "Ciliated Epithelial 2", '9' = "Blood Endothelial", '10' = "Pericyte", '11' = "Intermediate Epithelial", '12' = "T/NK Cell", '13' = "Epithelial/Fibroblast", '14' = "Macrophage", '15' = "Erythrocyte", '16' = "Luteal",'17' = "Mesothelial",'18' = "Lymph Endothelial/Epithelial") # Remove cluster due few data points and suspected doubletDistal_Named@active.ident <- factor(x = Distal_Named@active.ident, levels = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Pericyte','Blood Endothelial','Lymph Endothelial/Epithelial','Epithelial/Fibroblast','Stem-like Epithelial 1','Stem-like Epithelial 2','Intermediate Epithelial','Secretory Epithelial','Ciliated Epithelial 1','Ciliated Epithelial 2','T/NK Cell','Macrophage','Erythrocyte','Mesothelial','Luteal'))Distal_Named <- SetIdent(Distal_Named, value = Distal_Named@active.ident)Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
Muscle <- c('#E55451' , '#FFB7B2') # Reds
Endothelial <- c('#A0E6FF')  # Reds
FiboEpi <- "#FFE0B3" # Reddish Brown
Epi <-c('#6E3E6E','#8A2BE2','#CCCCFF','#DA70D6','#DF73FF','#604791') # Blues/Purples
Immune <- c( '#5A5E6B'  , '#B8C2CC' , '#FC86AA') # Yellowish Brown
Meso <- "#9EFFFF" # Pink
Lut <- "#9DCC00" # Greencolors <- c(Fibroblasts, Muscle, Endothelial, FiboEpi, Epi, Immune, Meso, Lut)pie(rep(1,length(colors)), col=colors) Distal_Named <- subset(Distal_Named, idents = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Pericyte','Blood Endothelial','Epithelial/Fibroblast','Stem-like Epithelial 1','Stem-like Epithelial 2','Intermediate Epithelial','Secretory Epithelial','Ciliated Epithelial 1','Ciliated Epithelial 2','T/NK Cell','Macrophage','Erythrocyte','Mesothelial','Luteal'))p1 <- DimPlot(Distal_Named,reduction='umap',cols=colors,pt.size = 0.5,label.size = 4,label.color = "black",repel = TRUE,label=F) +NoLegend() +labs(x="UMAP_1",y="UMAP_2")LabelClusters(p1, id="ident", color = "black", repel = T , size = 4, box.padding = .75)ggsave(filename = "FIG1b_all_distal_umap.pdf", plot = p1, width = 8, height = 12, dpi = 600)## Figure 1c: Distal Uterine Tube Features for Cell Type Identification ##features <- c("Dcn","Col1a1",           # Fibroblasts        "Acta2","Myh11","Myl9",   # Smooth Muscle"Pdgfrb","Mcam","Cspg4",  # Pericyte"Sele","Vwf","Tek",             # Blood Endothelial"Lyve1","Prox1","Icam1",          # Lymph Endothelial"Epcam","Krt8",           # Epithelial"Foxj1",                  # Ciliated Epithelial"Ovgp1",                  # Secretory Epithelial"Slc1a3","Pax8","Cd44","Itga6",         # Stem-like Epithelieal "Ptprc",                  # Immune                            "Cd8a","Cd4","Cd3e",      # T-Cell                            "Klrc1","Runx3",          # T/NK Cell"Klrd1",                  # NK Cell"Aif1","Cd68","Csf1r","Itgax", # Macrophage"Hbb-bs", "Hba-a1",       # Erythrocytes"Fras1","Rspo1","Lrrn4","Msln", # Mesothelial"Cyp11a1","Bcat1","Fkbp5","Spp1","Prlr") # Luteal Cellsall_dp <- DotPlot(object = Distal_Named,                    # Seurat objectassay = 'RNA',                        # Name of assay to use.  Default is the active assayfeatures = features,                  # List of features (select one from above or create a new one)# Colors to be used in the gradientcol.min = 0,                       # Minimum scaled average expression threshold (everything smaller will be set to this)col.max = 2.5,                        # Maximum scaled average expression threshold (everything larger will be set to this)dot.min = 0,                          # The fraction of cells at which to draw the smallest dot (default is 0)dot.scale = 9,                        # Scale the size of the pointsgroup.by = NULL,              # How the cells are going to be groupedsplit.by = NULL,                      # Whether to split the data (if you fo this make sure you have a different color for each variable)scale = TRUE,                         # Whether the data is scaledscale.by = "radius",                  # Scale the size of the points by 'size' or 'radius'scale.min = NA,                       # Set lower limit for scalingscale.max = NA                        # Set upper limit for scaling
)+    labs(x = NULL, y = NULL)+scale_color_viridis_c(option="F",begin=.4,end=0.9, direction = -1)+geom_point(aes(size=pct.exp), shape = 21, colour="black", stroke=0.6)+#theme_linedraw()+guides(x =  guide_axis(angle = 90))+ theme(axis.text.x = element_text(size = 14 , face = "italic"))+theme(axis.text.y = element_text(size = 14))+scale_y_discrete(limits = rev(levels(Distal_Named)))+theme(legend.title = element_text(size = 14))ggsave(filename = "FIG1c_all_distal_dotplot.pdf", plot = all_dp, width = 18, height = 10, dpi = 600)

图2

#### Figure 2: Characterization of distal epithelial cell states ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)library(monocle3)#### Load Distal Epithelial Dataset ####Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))EpiSecStemMarkers <- FindMarkers(Epi_Named, ident.1 = "Spdef+ Secretory",ident.2 = "Slc1a3+ Stem/Progenitor")
write.csv(EpiSecStemMarkers, file = "20240319_Staining_Markers2.csv")#### Figure 2a: Epithelial Cells of the Distal Uterine Tube ####epi_umap <- DimPlot(object = Epi_Named,                # Seurat object  reduction = 'umap',                 # Axes for the plot (UMAP, PCA, etc.) repel = TRUE,                       # Whether to repel the cluster labelslabel = FALSE,                       # Whether to have cluster labels cols = c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6"), #9pt.size = 0.6,                      # Size of each dot is (0.1 is the smallest)label.size = 0.5) +                   # Font size for labels    # You can add any ggplot2 1customizations herelabs(title = 'Colored by Cluster')+        # Plot titleNoLegend() +labs(x="UMAP_1",y="UMAP_2")ggsave(filename = "Fig2a_epi_umap.pdf", plot = epi_umap, width = 15, height = 12, dpi = 600)#### Figure 2c: Distal Uterine Tube Features for Epithelial Cell State Identification ####distal_features <- c("Krt8","Epcam","Slc1a3","Cd44","Sox9","Ovgp1","Sox17","Pax8", "Egr1","Itga6", "Bmpr1b","Rhoj", "Klf6","Msln","Cebpd","Dpp6", "Sec14l3", "Fam161a","Prom1", "Ly6a", "Kctd8", "Adam8","Dcn", "Col1a1", "Col1a2", "Timp3", "Pdgfra","Lgals1","Upk1a", "Thrsp","Spdef","Lcn2","Selenop", "Gstm2","Foxj1","Fam183b","Rgs22","Dnali1", "Mt1" , "Dynlrb2","Cdh1")epi_dp <- DotPlot(object = Epi_Named,                    # Seurat objectassay = 'RNA',                        # Name of assay to use.  Default is the active assayfeatures = distal_features,                  # List of features (select one from above or create a new one)# Colors to be used in the gradientcol.min = 0,                       # Minimum scaled average expression threshold (everything smaller will be set to this)col.max = 2.5,                        # Maximum scaled average expression threshold (everything larger will be set to this)dot.min = 0,                          # The fraction of cells at which to draw the smallest dot (default is 0)dot.scale = 4,                        # Scale the size of the pointsgroup.by = NULL,              # How the cells are going to be groupedsplit.by = NULL,                      # Whether to split the data (if you fo this make sure you have a different color for each variable)scale = TRUE,                         # Whether the data is scaledscale.by = "radius",                  # Scale the size of the points by 'size' or 'radius'scale.min = NA,                       # Set lower limit for scalingscale.max = NA )+                       # Set upper limit for scalinglabs(x = NULL,                              # x-axis labely = NULL)+scale_color_viridis_c(option="F",begin=.4,end=0.9, direction = -1)+geom_point(aes(size=pct.exp), shape = 21, colour="black", stroke=0.6)+#theme_linedraw()+guides(x =  guide_axis(angle = 90))+theme(axis.text.x = element_text(size = 8 , face = "italic"))+theme(axis.text.y = element_text(size = 9))+theme(legend.title = element_text(size = 9))+theme(legend.text = element_text(size = 8))+ scale_y_discrete(limits = c("Ciliated 2","Ciliated 1","Selenop+/Gstm2high Secretory","Spdef+ Secretory","Fibroblast-like","Pax8low/Prom1+ Cilia-forming", "Slc1a3med/Sox9+ Cilia-forming","Cebpdhigh/Foxj1- Progenitor","Slc1a3+ Stem/Progenitor"))ggsave(filename = "Fig2b_epi_dot_plot.pdf", plot = epi_dp, width = 8.3, height = 4.0625, dpi = 600)#### Load Distal Epithelial Pseudotime Dataset ####Distal_PHATE <- readRDS(file = "../dataset/Distal_Epi_PHATE.rds" , refhook = NULL)Beeg_PHATE <- Distal_PHATEBeeg_PHATE@reductions[["phate"]]@cell.embeddings <- Distal_PHATE@reductions[["phate"]]@cell.embeddings*100cds <- readRDS(file = "../dataset/Distal_Epi_PHATE_Monocle3.rds" , refhook = NULL)#### Figure 2b: PHATE embedding for differentiation trajectory of distal epithelial cells ####phate_dif <- DimPlot(Beeg_PHATE , reduction = "phate", cols = c("#B20224", "#35EFEF", "#00A1C6", "#A374B5", "#9000C6", "#EA68E1", "#59D1AF", "#2188F7", "#F28D86"),pt.size = 0.7,shuffle = TRUE,seed = 0,label = FALSE)+  labs(title = 'Colored by Cluster')+        # Plot titleNoLegend() +labs(x="UMAP_1",y="UMAP_2")ggsave(filename = "Fig3a_epi_phate.pdf", plot = phate_dif, width = 15, height = 12, dpi = 600)#### Figure 2d: PHATE and Monocle3 differentiation trajectory path ####pseudtotime <- plot_cells(cds, color_cells_by = "pseudotime",label_cell_groups=FALSE,label_leaves=FALSE,label_branch_points=FALSE,graph_label_size=0,cell_size = .01,cell_stroke = 1)+theme(axis.title.x = element_blank())+theme(axis.title.y = element_blank())+theme(axis.line.x = element_blank())+theme(axis.line.y = element_blank())+theme(axis.ticks.x = element_blank())+theme(axis.ticks.y = element_blank())+theme(axis.text.x = element_blank())+theme(axis.text.y = element_blank())+theme(legend.text = element_text(size = 12))ggsave(filename = "Fig3b_epi_pseudtotime.pdf", plot = pseudtotime, width = 18, height = 12, dpi = 600)#### Figure 2e: Slc1a3 PHATE Feature Plot ####Slc1a3_PHATE <- FeaturePlot(Beeg_PHATE, features = c("Slc1a3"), reduction = "phate", pt.size = 0.5)+scale_color_viridis_c(option="F",begin=.4,end=0.99, direction = -1)+theme(plot.title = element_text(size = 32,face = "bold.italic"))+theme(axis.title.x = element_blank())+theme(axis.title.y = element_blank())+theme(axis.line.x = element_blank())+theme(axis.line.y = element_blank())+theme(axis.ticks.x = element_blank())+theme(axis.ticks.y = element_blank())+theme(axis.text.x = element_blank())+theme(axis.text.y = element_blank())+theme(legend.text = element_text(size = 12))ggsave(filename = "Fig3c_SLC1A3_PHATE.pdf", plot = Slc1a3_PHATE, width = 18, height = 9, dpi = 600)#### Figure 2f: Pax8 PHATE Feature Plot ####Pax8_PHATE <- FeaturePlot(Beeg_PHATE, features = c("Pax8"), reduction = "phate", pt.size = 0.5)+scale_color_viridis_c(option="F",begin=.4,end=0.99, direction = -1)+theme(plot.title = element_text(size = 32,face = "bold.italic"))+theme(axis.title.x = element_blank())+theme(axis.title.y = element_blank())+theme(axis.line.x = element_blank())+theme(axis.line.y = element_blank())+theme(axis.ticks.x = element_blank())+theme(axis.ticks.y = element_blank())+theme(axis.text.x = element_blank())+theme(axis.text.y = element_blank())+theme(legend.text = element_text(size = 12))ggsave(filename = "Fig3d_PAX8_PHATE.pdf", plot = Pax8_PHATE, width = 18, height = 9, dpi = 600)

图6

#### Figure 6: Identification of cancer-prone cell states ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)library(monocle3)
library(ComplexHeatmap)
library(ggExtra)
library(gridExtra)
library(egg)library(scales)#### Load and Prepare Distal Epithelial and Epithelial Pseudotime Datasets ####Distal_PHATE <- readRDS(file = "../dataset/Distal_Epi_PHATE.rds" , refhook = NULL)cds <- readRDS(file = "../dataset/Distal_Epi_PHATE_Monocle3.rds" , refhook = NULL)Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))## Calculate Pseudotime Values ##pseudo <- pseudotime(cds)Distal_PHATE@meta.data$Pseudotime <- pseudo # Add to Seurat Metadata## Subset Seurat Object ##color_cells <- DimPlot(Distal_PHATE , reduction = "phate", cols = c("#B20224", #1"#35EFEF", #2"#00A1C6", #3"#A374B5", #4"#9000C6", #5"#EA68E1", #6"lightgrey", #7"#2188F7", #8"#F28D86"),pt.size = 0.7,shuffle = TRUE,seed = 0,label = FALSE)## Psuedotime and Lineage Assignment ##cellID <- rownames(Distal_PHATE@reductions$phate@cell.embeddings)
phate_embeddings <- Distal_PHATE@reductions$phate@cell.embeddings
pseudotime_vals <- Distal_PHATE@meta.data$Pseudotimecombined_data <- data.frame(cellID, phate_embeddings, pseudotime_vals)# Calculate the Average PHATE_1 Value for Pseudotime Points = 0 #
avg_phate_1 <- mean(phate_embeddings[pseudotime_vals == 0, 1])# Pseudotime Values lower than avge PHATE_1 Embedding will be Negative to split lineages
combined_data$Split_Pseudo <- ifelse(phate_embeddings[, 1] < avg_phate_1, -pseudotime_vals, pseudotime_vals)# Define Lineage #
combined_data$lineage <- ifelse(combined_data$PHATE_1 < avg_phate_1, "Secretory",ifelse(combined_data$PHATE_1 > avg_phate_1, "Ciliogenic", "Progenitor"))Distal_PHATE$Pseudotime_Adj <- combined_data$Split_Pseudo
Distal_PHATE$Lineage <- combined_data$lineage# Subset #Pseudotime_Lineage <- subset(Distal_PHATE, idents = c("Secretory 1","Secretory 2","Msln+ Progenitor","Slc1a3+/Sox9+ Cilia-forming","Pax8+/Prom1+ Cilia-forming","Progenitor","Ciliated 1","Ciliated 2"))## Set Bins ##bins <- cut_number(Pseudotime_Lineage@meta.data$Pseudotime_Adj , 40) # Evenly distribute bins Pseudotime_Lineage@meta.data$Bin <- bins # Metadata for Bins## Set Idents to PSeudoime Bin ##time_ident <- SetIdent(Pseudotime_Lineage, value = Pseudotime_Lineage@meta.data$Bin)av.exp <- AverageExpression(time_ident, return.seurat = T)$RNA # Calculate Avg log normalized expression# Calculates Average Expression for Each Bin #
# if you set return.seurat=T, NormalizeData is called which by default performs log-normalization #
# Reported as avg log normalized expression ##### Figure 6c: PHATE embedding for differentiation trajectory of distal epithelial cells ##### Create the stacked barplot
rainbow20 <- c('#FF0000','#FF6000','#FF8000','#FFA000','#FFC000','#FFE000','#FFFF00','#E0FF00','#C0FF00','#A0FF00','#00FF00','#00FFA0','#00F0FF','#00A0FF','#0020FF','#4000FF','#8000FF','#A000FF','#C000FF','#E000FF')rainbow_pseudo <- DimPlot(Pseudotime_Lineage , reduction = "phate", cols = c(rev(rainbow20),rainbow20),pt.size = 1.2,shuffle = TRUE,seed = 0,label = FALSE,group.by = "Bin")+    NoLegend()ggsave(filename = "rainbow_pseudo.pdf", plot = rainbow_pseudo, width = 20, height = 10, dpi = 600)#### Figure 6d: PHATE embedding for differentiation trajectory of distal epithelial cells ###### Pseudotime Scale Bar ##list <- 1:40
colors = c(rev(rainbow20),rainbow20)
df <- data.frame(data = list, color = colors)pseudo_bar <- ggplot(df, aes(x = 1:40, y = 1, fill = color)) + geom_bar(stat = "identity",position = "fill", color = "black", size = 0, width = 1) +scale_fill_identity() +theme_void()+ theme(axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())ggsave(filename = "pseudo_bar.pdf", plot = pseudo_bar, width = 0.98, height = 0.19, dpi = 600)## Plot gene list across pseudotime bin ##features <- c('Upk1a', "Spdef", "Ovgp1", "Gstm2", "Selenop", "Msln", "Slc1a3","Itga6", "Pax8",'Ecrg4', 'Sox5', 'Pde4b', 'Lcn2','Klf6','Trp53' , 'Trp73','Krt5','Foxa2','Prom1','Clstn2','Spef2','Dnah12','Foxj1', 'Fam166c' , 'Cfap126','Fam183b')# Create Bin List and expression of features #bin_list <- unique(Pseudotime_Lineage@meta.data$Bin) plot_info <- as.data.frame(av.exp[features, ]) # Call Avg Expression for featuresz_score <- transform(plot_info, SD=apply(plot_info,1, mean, na.rm = TRUE))
z_score <- transform(z_score, MEAN=apply(plot_info,1, sd, na.rm = TRUE))z_score1 <- (plot_info-z_score$MEAN)/z_score$SDplot_info$y <- rownames(plot_info) # y values as features
z_score1$y <- rownames(plot_info)plot_info <- gather(data = plot_info, x, expression, bin_list) #set plot
z_score1 <- gather(data = z_score1, x, z_score, bin_list) #set plot# Create Cell Clusters DF #Labeled_Pseudotime_Lineage <- RenameIdents(Pseudotime_Lineage, 'Secretory 1' = "Spdef+ Secretory", 'Progenitor' = "Slc1a3+ Stem/Progenitor", 'Msln+ Progenitor' = "Cebpdhigh/Foxj1- Progenitor",'Ciliated 1' = "Ciliated 1", 'Ciliated 2' = "Ciliated 2", 'Pax8+/Prom1+ Cilia-forming' = "Pax8low/Prom1+ Cilia-forming", 'Fibroblast-like' = "Fibroblast-like", #removed'Slc1a3+/Sox9+ Cilia-forming' = "Slc1a3med/Sox9+ Cilia-forming",'Secretory 2' = "Selenop+/Gstm2high Secretory")cluster_table <- table(Labeled_Pseudotime_Lineage@active.ident, Labeled_Pseudotime_Lineage@meta.data$Bin)clusters <- data.frame(cluster_table)clusters <- clusters %>% group_by(Var2) %>%mutate(Perc = Freq / sum(Freq))# Create Pseudotime DF #pseudotime_table <- table(seq(1, length(bin_list), 1), unique(Labeled_Pseudotime_Lineage@meta.data$Bin),seq(1, length(bin_list), 1))pseudotime_bins <- data.frame(pseudotime_table)  # calculate max and min z-scores
max_z <- max(z_score1$z_score, na.rm = TRUE)
min_z <- min(z_score1$z_score, na.rm = TRUE)# set color for outliers
outlier_color <- ifelse(z_score1$z_score > max_z | z_score1$z_score < min_z, ifelse(z_score1$z_score > 0, "#AD1F24", "#51A6DC"), "#e2e2e2")## Plot Gene Expression ### Set different na.value options for positive and negative values
na_color_pos <- "#AD1F24" # color for positive NA values
na_color_neg <- "#51A6DC" # color for negative NA valuescustom_bin_names <- c(paste0("S", 20:1), paste0("C", 1:20))figure <- ggplot(z_score1, aes(x, y, fill = z_score)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradientn(colors=c("#1984c5", "#e2e2e2", "#c23728"), name = "Average Expression \nZ-Score", limits = c(-3,3), na.value = ifelse(is.na(z_score1) & z_score1 > 0, na_color_pos, ifelse(is.na(z_score1) & z_score1 < 0, na_color_neg, "grey50")),oob = scales::squish)+scale_x_discrete(limits= sort(bin_list) , labels= custom_bin_names)+scale_y_discrete(limits= rev(features))+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 0, hjust = 0.5, size = 10, face = "bold"),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold.italic"))+theme(plot.title = element_blank(),plot.margin=unit(c(-0.5,1,1,1), "cm"))## Plot Cluster Percentage ##`Spdef+ Secretory` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Spdef+ Secretory")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(1,1,1,1), "cm"))`Selenop+/Gstm2high Secretory` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Selenop+/Gstm2high Secretory")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Cebpdhigh/Foxj1- Progenitor` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Cebpdhigh/Foxj1- Progenitor")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Slc1a3+ Stem/Progenitor` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Slc1a3+ Stem/Progenitor")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Slc1a3med/Sox9+ Cilia-forming` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Slc1a3med/Sox9+ Cilia-forming")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Pax8low/Prom1+ Cilia-forming` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Pax8low/Prom1+ Cilia-forming")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Ciliated 1` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Ciliated 1")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Ciliated 2` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Ciliated 2")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))## Plot Pseudotime Color ##list <- 1:40
colors = c(rev(rainbow20),rainbow20)
df <- data.frame(data = list, color = colors)binning <- ggplot(df, aes(x = 1:40, y = 1, fill = color)) + geom_bar(stat = "identity",position = "fill", color = "black", size = 1, width = 1) +scale_fill_identity() +theme_void()+ theme(axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Pseudotime Bin ")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust =1, vjust = .75, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))### Combine Plots ###psuedotime_lineage <- ggarrange(`Spdef+ Secretory`,`Selenop+/Gstm2high Secretory`,`Cebpdhigh/Foxj1- Progenitor`,`Slc1a3+ Stem/Progenitor`,`Slc1a3med/Sox9+ Cilia-forming`,`Pax8low/Prom1+ Cilia-forming`,`Ciliated 1`,`Ciliated 2`,`binning`,figure , ncol=1,heights = c(2, 2, 2, 2, 2, 2, 2, 2, 2, (2*length(features)),widths = c(3)),padding = unit(0.01))ggsave(filename = "FIG6d_psuedotime_lineage.pdf", plot = psuedotime_lineage, width = 18, height = 9, dpi = 600)#### Figure 6e: Stacked violin plots for cancer-prone cell states ####### Stacked Violin Plot Function ####https://divingintogeneticsandgenomics.rbind.io/post/stacked-violin-plot-for-visualizing-single-cell-data-in-seurat/## remove the x-axis text and tick
## plot.margin to adjust the white space between each plot.
## ... pass any arguments to VlnPlot in Seurat
modify_vlnplot <- function(obj, feature, pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {p<- VlnPlot(obj, features = feature, pt.size = pt.size, ... )  + xlab("") + ylab(feature) + ggtitle("") + theme(legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = rel(1), angle = 0, face = "bold.italic"), axis.text.y = element_text(size = rel(1)), plot.margin = plot.margin ) return(p)
}## extract the max value of the y axis
extract_max<- function(p){ymax<- max(ggplot_build(p)$layout$panel_scales_y[[1]]$range$range)return(ceiling(ymax))
}## main function
StackedVlnPlot<- function(obj, features,pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {plot_list<- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))# Add back x-axis title to bottom plot. patchwork is going to support this?plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +theme(axis.text.x=element_text(angle = 60, hjust=1, vjust=0.95), axis.ticks.x = element_line())# change the y-axis tick to only max value ymaxs<- purrr::map_dbl(plot_list, extract_max)plot_list<- purrr::map2(plot_list, ymaxs, function(x,y) x + scale_y_continuous(breaks = c(y)) + expand_limits(y = y))p<- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)return(p)
}features<- c("Slc1a3", "Pax8" , "Trp73" , "Prom1" )features<- c("Pax8", "Ovgp1" ,  "Lcn2" , "Upk1a" , "Spdef" ,"Thrsp" )colors <- c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4#"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6")
No_Fibro <- subset(x = Epi_Named, idents =  c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", #"Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2"))stack_vln <- StackedVlnPlot(obj = No_Fibro, features = features, slot = "data",pt.size = 0,cols = c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4#"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6"))+ #9theme(plot.title = element_text(size = 32, face = "bold.italic"))+scale_x_discrete(limits = c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", #"Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2"))+theme(axis.text.x = element_text(size = 16, angle = 60))+theme(axis.text.y = element_text(size = 14))+theme(axis.title.y.left = element_text(size = 16))ggsave(filename = "FIG6e_stacked_vln_noFibro.pdf", plot = stack_vln, width = 18, height = 12, dpi = 600)#### Figure 6f: Krt5 expression within epithelial cell states ####ggsave(filename = "Krt5_dp_others.pdf", plot = Krt5_dp_others, width = 1.89*8, height = 3.06*8, dpi = 600)Krt5_dp <- DotPlot(object = No_Fibro,                    # Seurat objectassay = 'RNA',                        # Name of assay to use.  Default is the active assayfeatures = 'Krt5',                  # List of features (select one from above or create a new one)# Colors to be used in the gradientcol.min = 0,                       # Minimum scaled average expression threshold (everything smaller will be set to this)col.max = 2.5,                        # Maximum scaled average expression threshold (everything larger will be set to this)dot.min = 0,                          # The fraction of cells at which to draw the smallest dot (default is 0)dot.scale = 24,                        # Scale the size of the pointsgroup.by = NULL,              # How the cells are going to be groupedsplit.by = NULL,                      # Whether to split the data (if you fo this make sure you have a different color for each variable)scale = TRUE,                         # Whether the data is scaledscale.by = "radius",                  # Scale the size of the points by 'size' or 'radius'scale.min = NA,                       # Set lower limit for scalingscale.max = NA )+                       # Set upper limit for scalinglabs(x = NULL,                              # x-axis labely = NULL)+scale_color_viridis_c(option="F",begin=.4,end=0.9, direction = -1)+geom_point(aes(size=pct.exp), shape = 21, colour="black", stroke=0.7)+theme_linedraw(base_line_size = 5)+guides(x =  guide_axis(angle = 90))+theme(axis.text.x = element_text(size = 32 , face = "italic"))+theme(axis.text.y = element_text(size = 32))+theme(legend.title = element_text(size = 12))+  scale_y_discrete(limits = c("Ciliated 2","Ciliated 1","Selenop+/Gstm2high Secretory","Spdef+ Secretory",#"Fibroblast-like","Pax8low/Prom1+ Cilia-forming", "Slc1a3med/Sox9+ Cilia-forming","Cebpdhigh/Foxj1- Progenitor","Slc1a3+ Stem/Progenitor"))ggsave(filename = "Krt5_dp_noFibro.pdf", plot = Krt5_dp, width = 1.89*8, height = 3.06*8, dpi = 600)

附图2

#### Figure Supp 2: Characterization of distal epithelial cell states ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)#### Unprocessed and Processed Distal Dataset ####All.merged <- readRDS(file = "../dataset/Unfiltered_Mouse_Distal.rds", refhook = NULL) # Prior to Quality ControlDistal <- readRDS(file = "../dataset/Distal_Filtered_Cells.rds" , refhook =  NULL) # After to Quality ControlDistal_Named <- RenameIdents(Distal, '0' = "Fibroblast 1", '1' = "Fibroblast 2", '2' = "Secretory Epithelial",'3' = "Smooth Muscle", '4' = "Ciliated Epithelial 1", '5' = "Fibroblast 3", '6' = "Stem-like Epithelial 1",'7' = "Stem-like Epithelial 2",'8' = "Ciliated Epithelial 2", '9' = "Blood Endothelial", '10' = "Pericyte", '11' = "Intermediate Epithelial", '12' = "T/NK Cell", '13' = "Epithelial/Fibroblast", '14' = "Macrophage", '15' = "Erythrocyte", '16' = "Luteal",'17' = "Mesothelial",'18' = "Lymph Endothelial/Epithelial") # Remove cluster due few data points and suspected doubletDistal_Named@active.ident <- factor(x = Distal_Named@active.ident, levels = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Pericyte','Blood Endothelial','Lymph Endothelial/Epithelial','Epithelial/Fibroblast','Stem-like Epithelial 1','Stem-like Epithelial 2','Intermediate Epithelial','Secretory Epithelial','Ciliated Epithelial 1','Ciliated Epithelial 2','T/NK Cell','Macrophage','Erythrocyte','Mesothelial','Luteal'))Distal_Named <- subset(Distal_Named, idents = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Pericyte','Blood Endothelial','Epithelial/Fibroblast','Stem-like Epithelial 1','Stem-like Epithelial 2','Intermediate Epithelial','Secretory Epithelial','Ciliated Epithelial 1','Ciliated Epithelial 2','T/NK Cell','Macrophage','Erythrocyte','Mesothelial','Luteal'))Distal_Named <- SetIdent(Distal_Named, value = Distal_Named@active.ident)Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))#### Figure Supp 2a: Unfilitered % MT Genes ####unfiltered_MT <- VlnPlot(All.merged, features = c("percent.mt"), group.by = 'Sample', pt.size = 0,cols = natparks.pals(name="Arches2",n=3))+theme(legend.position = 'none')+theme(axis.text.x = element_text(size = 16))+     # Change X-Axis Text Sizetheme(axis.text.y = element_text(size = 16))+     # Change Y-Axis Text Sizetheme(axis.title.y = element_text(size = 18))+    # Change Y-Axis Title Text Sizetheme(plot.title = element_text(size = 32,face = "bold.italic"))+theme(axis.title.x = element_blank())ggsave(filename = "FIGs2a_unfiltered_MT.pdf", plot = unfiltered_MT, width = 12, height = 9, dpi = 600)#exprs <- as.data.frame(FetchData(object = All.merged, vars = c('nCount_RNA' , "Sample")))#df_new <- filter(exprs, Sample == 'mD1')#mean(df_new$nCount_RNA)
#sd(df_new$nCount_RNA)# Remove the original 'Value' column if needed
df_new <- df_new %>%select(-Value)x <- spread(exprs, Sample, percent.mt )
mean(exprs)
#### Figure Supp 2b: Unfilitered nFeature RNA #### unfiltered_nFeature <- VlnPlot(All.merged, features = c("nFeature_RNA"), group.by = 'Sample', pt.size = 0,cols = natparks.pals(name="Arches2",n=3))+theme(legend.position = 'none')+theme(axis.text.x = element_text(size = 16))+theme(axis.text.y = element_text(size = 16))+theme(axis.title.y = element_text(size = 18))+theme(plot.title = element_text(size = 32,face = "bold.italic"))+theme(axis.title.x = element_blank()) # Change object to visualize other samplesggsave(filename = "FIGs2b_unfiltered_nFeature.pdf", plot = unfiltered_nFeature, width = 12, height = 9, dpi = 600)#### Figure Supp 2c: Unfilitered nCount RNA #### unfiltered_nCount <- VlnPlot(All.merged, features = c("nCount_RNA"), group.by = 'Sample', pt.size = 0,cols = natparks.pals(name="Arches2",n=3))+theme(legend.position = 'none')+theme(axis.text.x = element_text(size = 16))+theme(axis.text.y = element_text(size = 16))+theme(axis.title.y = element_text(size = 18))+theme(plot.title = element_text(size = 32,face = "bold.italic"))+theme(axis.title.x = element_blank()) # Change object to visualize other samplesggsave(filename = "FIGs2c_unfiltered_nCount.pdf", plot = unfiltered_nCount, width = 12, height = 9, dpi = 600)#### Figure Supp 2d: Doublets All Cells #### All_Doublet <- DimPlot(object = Distal, reduction = 'umap', group.by = "Doublet",cols = c( "#ffb6c1", "#380b11"),repel = TRUE, label = F, pt.size = 1.2, order = c("Doublet","Singlet"),label.size = 5) +labs(x="UMAP_1",y="UMAP_2")ggsave(filename = "FIGs2d1_All_Doublet_umap.pdf", plot = All_Doublet, width = 22, height = 17, dpi = 600)## Stacked Bar Doublets ##table <- table(Distal_Named@active.ident ,Distal_Named@meta.data$Doublet)    # Create a table of countsdf <- data.frame(table) doublet <- ggplot(data = df,                # Dataset to use for plot.  Needs to be a data.frame  aes(x = Var1,              # Variable to plot on the x-axisy = Freq,              # Variable to plot on the y-axisfill = factor(Var2,    # Variable to fill the barslevels = c("Doublet","Singlet")))) + # Order of the stacked barstheme_classic() +               # ggplot2 theme# Bar plotgeom_bar(position = 'fill',     # Position of bars.  Fill means the bars are stacked.stat = "identity",     # Height of bars represent values in the datasize = 1) +            # Size of bars# Color schemescale_fill_manual("Doublet", limits = c("Doublet","Singlet"),values = c('#8B0000','#808080')) +# Add plot labelslabs(x = NULL,                     # x-axis labely = "Fraction of Cells") +    # y-axis labeltheme(text = element_text(size = 15),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 60, hjust = 1, size = 11),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1))+       # Text color and horizontal adjustment on y-axisscale_x_discrete(limits = (c('Intermediate Epithelial','Epithelial/Fibroblast','Stem-like Epithelial 1','Ciliated Epithelial 1','Erythrocyte','Smooth Muscle','Stem-like Epithelial 2','Mesothelial','Blood Endothelial','Pericyte','Fibroblast 2','Secretory Epithelial','Fibroblast 1','Ciliated Epithelial 2','Fibroblast 3','T/NK Cell','Macrophage','Luteal')))+coord_flip()ggsave(filename = "FIGs2d2_doublet_quant.pdf", plot = doublet, width = 10, height = 16, dpi = 600)#### Figure Supp 2e: Sample Distribution for All Cells #### table <- table(Distal_Named@active.ident ,Distal_Named@meta.data$Sample)    # Create a table of countsdf <- data.frame(table) table2 <- table(Epi_Named@meta.data$Sample)all_sample_dist <- ggplot(data = df,                # Dataset to use for plot.  Needs to be a data.frame  aes(x = Var1,              # Variable to plot on the x-axisy = Freq,              # Variable to plot on the y-axisfill = factor(Var2,    # Variable to fill the barslevels = c("mD1","mD2","mD4")))) + # Order of the stacked barstheme_classic() +               # ggplot2 theme# Bar plotgeom_bar(position = 'fill',     # Position of bars.  Fill means the bars are stacked.stat = "identity",     # Height of bars represent values in the datasize = 1) +            # Size of bars# Color schemescale_fill_manual("Location", limits = c("mD1","mD2","mD4"),values = c(natparks.pals(name="Arches2",n=3))) +# Add plot labelslabs(x = NULL,                     # x-axis labely = "Fraction of Cells") +    # y-axis labeltheme(text = element_text(size = 15),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 60, hjust = 1, size = 11),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1))               # Text color and horizontal adjustment on y-axisggsave(filename = "FIGs2f_all_sample_dist.pdf", plot = epi_sample_dist, width = 16, height = 12, dpi = 600)#### Figure Supp 2f: Doublets Epithelial Cells #### Epi_Doublet <- DimPlot(object = Epi_Named, reduction = 'umap', group.by = "Doublet",cols = c( "#ffb6c1", "#380b11"),repel = TRUE, label = F, pt.size = 1.2, order = c("Doublet","Singlet"),label.size = 5) +labs(x="UMAP_1",y="UMAP_2")ggsave(filename = "FIGs2e1_Epi_Doublet_umap.pdf", plot = All_Doublet, width = 22, height = 17, dpi = 600)## Stacked Bar Doublets ##table <- table(Epi_Named@active.ident ,Epi_Named@meta.data$Doublet)    # Create a table of countsdf <- data.frame(table) epi_doublet <- ggplot(data = df,                # Dataset to use for plot.  Needs to be a data.frame  aes(x = Var1,              # Variable to plot on the x-axisy = Freq,              # Variable to plot on the y-axisfill = factor(Var2,    # Variable to fill the barslevels = c("Doublet","Singlet")))) + # Order of the stacked barstheme_classic() +               # ggplot2 theme# Bar plotgeom_bar(position = 'fill',     # Position of bars.  Fill means the bars are stacked.stat = "identity",     # Height of bars represent values in the datasize = 1) +            # Size of bars# Color schemescale_fill_manual("Doublet", limits = c("Doublet","Singlet"),values = c('#8B0000','#808080')) +# Add plot labelslabs(x = NULL,                     # x-axis labely = "Fraction of Cells") +    # y-axis labeltheme(text = element_text(size = 15),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 60, hjust = 1, size = 11),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1))+       # Text color and horizontal adjustment on y-axisscale_x_discrete(limits = (c("Slc1a3med/Sox9+ Cilia-forming","Fibroblast-like","Pax8low/Prom1+ Cilia-forming","Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Ciliated 1","Ciliated 2", "Spdef+ Secretory","Selenop+/Gstm2high Secretory")))+coord_flip()ggsave(filename = "FIGs2e2_epi_doublet_quant.pdf", plot = epi_doublet, width = 10, height = 16, dpi = 600)#### Figure Supp 2g: Sample Distribution for Epithelial Cells #### table <- table(Epi_Named@active.ident ,Epi_Named@meta.data$Sample)    # Create a table of countsdf <- data.frame(table) table2 <- table(Epi_Named@meta.data$Samlpe)epi_sample_dist <- ggplot(data = df,                # Dataset to use for plot.  Needs to be a data.frame  aes(x = Var1,              # Variable to plot on the x-axisy = Freq,              # Variable to plot on the y-axisfill = factor(Var2,    # Variable to fill the barslevels = c("mD1","mD2","mD4")))) + # Order of the stacked barstheme_classic() +               # ggplot2 theme# Bar plotgeom_bar(position = 'fill',     # Position of bars.  Fill means the bars are stacked.stat = "identity",     # Height of bars represent values in the datasize = 1) +            # Size of bars# Color schemescale_fill_manual("Location", limits = c("mD1","mD2","mD4"),values = c(natparks.pals(name="Arches2",n=3))) +# Add plot labelslabs(x = NULL,                     # x-axis labely = "Fraction of Cells") +    # y-axis labeltheme(text = element_text(size = 15),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 60, hjust = 1, size = 11),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1))               # Text color and horizontal adjustment on y-axisggsave(filename = "FIGs2g_epi_sample_dist.pdf", plot = epi_sample_dist, width = 16, height = 12, dpi = 600)#### Figure Supp 2h: Distal Tile Mosaic ####library(treemap)dist_cell_types <- table(Idents(Distal_Named), Distal_Named$orig.ident)
dist_cell_type_df <- as.data.frame(dist_cell_types)## Colors ##Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
FiboEpi <- "#FFE0B3" # Reddish Brown
Muscle <- c('#E55451' , '#FFB7B2') # Reds
Endothelial <- c('#A0E6FF')  # Reds
Epi <-c('#6E3E6E','#8A2BE2','#604791','#CCCCFF','#DA70D6','#DF73FF') # Blues/Purples
Immune <- c( '#5A5E6B'  , '#B8C2CC' , '#FC86AA') # Yellowish Brown
Meso <- "#9EFFFF" # Pink
Lut <- "#9DCC00" # Greencolors <- c(Fibroblasts, FiboEpi, Muscle, Endothelial, Epi, Immune, Meso, Lut)## Tile Mosaic ##distal_treemap <- treemap(dist_cell_type_df, index = 'Var1', vSize= 'Freq', vColor = colors, palette = colors)ggsave(filename = "20240612_all_distal_tile.pdf", plot = distal_treemap, width = 12, height = 8, dpi = 600)#### Figure Supp 2i: Epi Markers All Distal Cells ####### Stacked Violin Plot Function ####https://divingintogeneticsandgenomics.rbind.io/post/stacked-violin-plot-for-visualizing-single-cell-data-in-seurat/## remove the x-axis text and tick
## plot.margin to adjust the white space between each plot.
## ... pass any arguments to VlnPlot in Seuratmodify_vlnplot <- function(obj, feature, pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {p<- VlnPlot(obj, features = feature, pt.size = pt.size, ... )  + xlab("") + ylab(feature) + ggtitle("") + theme(legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = rel(1), angle = 0, face = "bold.italic"), axis.text.y = element_text(size = rel(1)), plot.margin = plot.margin ) return(p)
}## extract the max value of the y axis
extract_max<- function(p){ymax<- max(ggplot_build(p)$layout$panel_scales_y[[1]]$range$range)return(ceiling(ymax))
}## main function
StackedVlnPlot<- function(obj, features,pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {plot_list<- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))# Add back x-axis title to bottom plot. patchwork is going to support this?plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +theme(axis.text.x=element_text(angle = 60, hjust=1, vjust=0.95), axis.ticks.x = element_line())# change the y-axis tick to only max value ymaxs<- purrr::map_dbl(plot_list, extract_max)# plot_list<- purrr::map2(plot_list, ymaxs, function(x,y) x + plot_list<- purrr::map2(plot_list, c(5,5,8,5), function(x,y) x +                           scale_y_continuous(breaks = c(y)) + expand_limits(y = y))p<- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)return(p)
}features<- c("Epcam", "Krt8" , "Ovgp1" , "Foxj1" )Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
Muscle <- c('#E55451' , '#FFB7B2') # Reds
Endothelial <- c('#A0E6FF')  # Reds
FiboEpi <- "#FFE0B3" # Reddish Brown
Epi <-c('#6E3E6E','#8A2BE2','#604791','#CCCCFF','#DA70D6','#DF73FF') # Blues/Purples
Immune <- c( '#5A5E6B'  , '#B8C2CC' , '#FC86AA') # Yellowish Brown
Meso <- "#9EFFFF" # Pink
Lut <- "#9DCC00" # Greencolors <- c(Fibroblasts, FiboEpi, Muscle, Endothelial, Epi, Immune, Meso, Lut)stack_vln <- StackedVlnPlot(obj = Distal_Named, features = features, slot = "data",pt.size = 0,cols = colors)+ #9theme(plot.title = element_text(size = 32, face = "bold.italic"))+scale_x_discrete(limits = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Epithelial/Fibroblast','Smooth Muscle','Pericyte','Blood Endothelial','Stem-like Epithelial 1','Stem-like Epithelial 2','Intermediate Epithelial','Secretory Epithelial','Ciliated Epithelial 1','Ciliated Epithelial 2','T/NK Cell','Macrophage','Erythrocyte','Mesothelial','Luteal'))+theme(axis.text.x = element_text(size = 16, angle = 60))+theme(axis.text.y = element_text(size = 14))+theme(axis.title.y.left = element_text(size = 16))ggsave(filename = "20240612_All_Distal_stacked_vln.pdf", plot = stack_vln, width = 18, height = 12, dpi = 600)#### Figure Supp 2j: Epi Markers Distal Epi Cells ####Epi_Sub <- subset(Epi_Named, idents = c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2"))colors <- c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6")stack_vln <- StackedVlnPlot(obj = Epi_Named, features = features, slot = "data",pt.size = 0,cols = c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6"))+ #9theme(plot.title = element_text(size = 32, face = "bold.italic"))+scale_x_discrete(limits = c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2"))+theme(axis.text.x = element_text(size = 16, angle = 60))+theme(axis.text.y = element_text(size = 14))+theme(axis.title.y.left = element_text(size = 16))ggsave(filename = "20240612_Distal_Epi_stacked_vln.pdf", plot = stack_vln, width = 18, height = 12, dpi = 600)

附图3

#### Figure Supp 3: Doublet detection of fibroblast and epithelial markers ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)#### Load Distal Epithelial Dataset ####Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))table(Epi_Named@active.ident)#### Plot Compilation ####feature_scatter <- FeatureScatter( Fibroblast, "Krt8","Col1a1",cols = c("#35EFEF", #1"#00A1C6", #2"#2188F7", #3"#EA68E1", #4"#59D1AF", #5"#B20224", #6"#F28D86", #7"#A374B5", #8"#9000C6"))x <- DotPlot(Epi_Named , features = c("Krt8" , "Col1a1"))
write.csv(x$data , "doublet_data.csv")Fibroblast <- subset(Epi_Named, idents = c("Fibroblast-like"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( Fibroblast, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Fibroblast_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Stem ##c("Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")Stem <- subset(Epi_Named, idents = c("Slc1a3+ Stem/Progenitor"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( Stem, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Slc1a3_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Prog ##c("Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")Prog <- subset(Epi_Named, idents = c("Cebpdhigh/Foxj1- Progenitor"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( Prog, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Cebpd_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Cilia-forming ##c("Pax8low/Prom1+ Cilia-forming","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")trans <- subset(Epi_Named, idents = c("Slc1a3med/Sox9+ Cilia-forming"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( trans, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_SlcCilia_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Pax8 Cilia-forming ##c("Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")cancer <- subset(Epi_Named, idents = c("Pax8low/Prom1+ Cilia-forming"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( cancer, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Prom1_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Spdef Secretory ##c("Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")sec1 <- subset(Epi_Named, idents = c("Spdef+ Secretory"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( sec1, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Spdef_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Selenop Secretory ##c("Ciliated 1","Ciliated 2")sec2 <- subset(Epi_Named, idents = c("Selenop+/Gstm2high Secretory"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( sec2, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Selenop_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Ciliated 1 ##c("Ciliated 2")cil1 <- subset(Epi_Named, idents = c("Ciliated 1"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( cil1, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Ciliated_1_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)## Ciliated 2 ##cil2 <- subset(Epi_Named, idents = c("Ciliated 2"))custom_labels <- function(x) {ifelse(x %% 1 == 0, as.character(x), "")
}feature_scatter <- FeatureScatter( cil2, "Krt8","Col1a1",cols = "black")+  # Scatter plotNoLegend()+labs(title = NULL)+theme(panel.grid.major = element_line(color = "grey", size = 0.5),panel.grid.minor = element_blank())+theme(axis.text.x = element_text(color = 'black', size = 12),axis.title.x = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'black'),axis.text.y = element_text(color = 'black', size = 12),axis.title.y = element_text(color = 'black', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'black'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5),labels = custom_labels) +  # Set breaks every 0.5 units on x-axisscale_y_continuous(limits = c(0,5) , breaks = seq(0, 10, by = 0.5),labels = custom_labels)   # Set breaks every 0.5 units on y-axisplot_data <- as.data.frame(feature_scatter$data)# Create density plot for x-axis
density_x <- ggplot(plot_data, aes(x = Krt8 , fill = 'black')) +geom_density() +theme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_blank(),axis.title.x = element_blank(),axis.ticks.x = element_blank(),axis.text.y = element_text(color = 'white', size = 12),axis.title.y = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.y = element_line(color = 'white'),plot.margin = unit(c(0, 0, 0, 0), "mm"))+scale_y_continuous(labels = function(y) sprintf("%.0f", y))+ NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,4) , breaks = seq(0, 10, by = 0.5))  # Set breaks every 0.5 units on x-axis# Create density plot for y-axis
density_y <- ggplot(plot_data, aes(x = Col1a1 , fill = 'black')) +geom_density() +coord_flip() +  # Flip axes for y-density plottheme(axis.line = element_line(color='white'),panel.background = element_blank()) +theme(axis.text.x = element_text(color = 'white', size = 12),axis.title.x = element_text(color = 'white', size = 14, face = "bold.italic"),axis.ticks.x = element_line(color = 'white'),axis.text.y = element_blank(),axis.title.y = element_blank(),axis.ticks.y = element_blank(),plot.margin = unit(c(0, 0, 0, 0), "mm"))+NoLegend()+scale_fill_manual(values = c("grey"))+scale_x_continuous(limits = c(0,5))+scale_y_continuous(limits = c(0,1) , breaks = seq(0, 10, by = 0.5))   # Set breaks every 0.5 units on y-axis# Arrange plotstop_row <- cowplot::plot_grid(density_x, NULL,nrow = 1, rel_widths = c(3, 1), rel_heights = c(0.5,0.5))bottom_row <- cowplot::plot_grid(feature_scatter,density_y,nrow = 1, rel_widths = c(3, 1), rel_heights = c(3,3))combined_plot <- cowplot::plot_grid(top_row , bottom_row,nrow = 2 , rel_widths = c(3, 1), rel_heights = c(0.25,3))ggsave(filename = "20240221_Ciliated_2_Doublet.pdf", plot = combined_plot, width = 12, height = 12, dpi = 600)

附图4

#### Figure Supp 4: Census of cell types of the mouse uterine tube ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)#### Proximal Datasets ####Proximal <- readRDS( file = "../dataset/Proximal_Filtered_Cells.rds" , refhook =  NULL)Proximal_Named <- RenameIdents(Proximal, '0' = "Fibroblast 1", '1' = "Stem-like Epithelial", '2' = "Fibroblast 2",'3' = "Fibroblast 3", '4' = "Immune", '5' = "Secretory Epithelial", '6' = "Endothelial",'7' = "Ciliated Epithelial",'8' = "Mesothelial", '9' = "Smooth Muscle")Proximal_Named@active.ident <- factor(x = Proximal_Named@active.ident, levels = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Endothelial','Stem-like Epithelial','Secretory Epithelial','Ciliated Epithelial','Immune','Mesothelial'))Proximal_Named <- SetIdent(Proximal_Named, value = Proximal_Named@active.ident)Epi_Filter <- readRDS(file = "../dataset/Proximal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Dbi+/Spdefhigh Secretory", '1' = "Bmpr1b+ Progenitor", '2' = "Wfdc2+ Secretory",'3' = "Ciliated", '4' = "Sox17high Secretory", '5' = "Kcne3+ Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c("Ciliated","Dbi+/Spdefhigh Secretory","Kcne3+ Secretory","Sox17high Secretory","Wfdc2+ Secretory","Bmpr1b+ Progenitor"))#### Figure Supp 4a: Proximal All Cell Types ####Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
Muscle <- c('#E55451') # Reds
Endothelial <- c('#A0E6FF')  # Blues
Epi <-c('#6E3E6E','#CCCCFF','#DF73FF') # Purples
Immune <- c( '#5A5E6B' ) # Grey
Meso <- "#1F51FF" # Neon BLuecolors <- c(Fibroblasts, Muscle, Endothelial, Epi, Immune, Meso)p1 <- DimPlot(Proximal_Named,reduction='umap',cols=colors,pt.size = 1.4,label.size = 4,label.color = "black",repel = TRUE,label=F) +NoLegend() +labs(x="UMAP_1",y="UMAP_2")ggsave(filename = "FIGs3a_all_proximal_umap.pdf", plot = p1, width = 15, height = 12, dpi = 600)#### Figure Supp 4b: Proximal Tile Mosaic  ####library(treemap)Proximal_Named <- RenameIdents(Proximal, '0' = "Fibroblast 1", '1' = "Stem-like Epithelial", '2' = "Fibroblast 2",'3' = "Fibroblast 3", '4' = "Immune", '5' = "Secretory Epithelial", '6' = "Endothelial",'7' = "Ciliated Epithelial",'8' = "Mesothelial", '9' = "Smooth Muscle")Proximal_Named@active.ident <- factor(x = Proximal_Named@active.ident, levels = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Endothelial','Stem-like Epithelial','Secretory Epithelial','Ciliated Epithelial','Immune','Mesothelial'))Proximal_Named <- SetIdent(Proximal_Named, value = Proximal_Named@active.ident)Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
Muscle <- c('#E55451') # Reds
Endothelial <- c('#A0E6FF')  # Blues
Epi <-c('#6E3E6E','#CCCCFF','#DF73FF') # Purples
Immune <- c( '#5A5E6B' ) # Grey
Meso <- "#1F51FF" # Neon BLuecolors <- c(Fibroblasts, Muscle, Endothelial, Epi, Immune, Meso)prox_cell_types <- table(Idents(Proximal_Named), Proximal_Named$orig.ident)
prox_cell_type_df <- as.data.frame(prox_cell_types)## Tile Mosaic ##prox_treemap <- treemap(prox_cell_type_df, index = 'Var1', vSize= 'Freq', vColor = colors, palette = colors)ggsave(filename = "20240612_all_prox_tile.pdf", plot = prox_cell_type_df, width = 8, height = 12, dpi = 600)#### Figure Supp 4c: Epi Markers All Distal Cells ####### Stacked Violin Plot Function ####https://divingintogeneticsandgenomics.rbind.io/post/stacked-violin-plot-for-visualizing-single-cell-data-in-seurat/## remove the x-axis text and tick
## plot.margin to adjust the white space between each plot.
## ... pass any arguments to VlnPlot in Seuratmodify_vlnplot <- function(obj, feature, pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {p<- VlnPlot(obj, features = feature, pt.size = pt.size, ... )  + xlab("") + ylab(feature) + ggtitle("") + theme(legend.position = "none", axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = rel(1), angle = 0, face = "bold.italic"), axis.text.y = element_text(size = rel(1)), plot.margin = plot.margin ) return(p)
}## extract the max value of the y axis
extract_max<- function(p){ymax<- max(ggplot_build(p)$layout$panel_scales_y[[1]]$range$range)return(ceiling(ymax))
}## main function
StackedVlnPlot<- function(obj, features,pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {plot_list<- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))# Add back x-axis title to bottom plot. patchwork is going to support this?plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +theme(axis.text.x=element_text(angle = 60, hjust=1, vjust=0.95), axis.ticks.x = element_line())# change the y-axis tick to only max value ymaxs<- purrr::map_dbl(plot_list, extract_max)# plot_list<- purrr::map2(plot_list, ymaxs, function(x,y) x + plot_list<- purrr::map2(plot_list, c(5,5,8,5), function(x,y) x +                           scale_y_continuous(breaks = c(y)) + expand_limits(y = y))p<- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)return(p)
}features<- c("Epcam", "Krt8" , "Ovgp1" , "Foxj1" )Fibroblasts <- c('#FF9D00' , '#FFB653' , '#FFCB9A')   # Oranges
Muscle <- c('#E55451') # Reds
Endothelial <- c('#A0E6FF')  # Blues
Epi <-c('#6E3E6E','#CCCCFF','#DF73FF') # Purples
Immune <- c( '#5A5E6B' ) # Grey
Meso <- "#1F51FF" # Neon BLuecolors <- c(Fibroblasts, Muscle, Endothelial, Epi, Immune, Meso)stack_vln <- StackedVlnPlot(obj = Proximal_Named, features = features, slot = "data",pt.size = 0,cols = colors)+ #9theme(plot.title = element_text(size = 32, face = "bold.italic"))+scale_x_discrete(limits = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Endothelial','Stem-like Epithelial','Secretory Epithelial','Ciliated Epithelial','Immune','Mesothelial'))+theme(axis.text.x = element_text(size = 16, angle = 60))+theme(axis.text.y = element_text(size = 14))+theme(axis.title.y.left = element_text(size = 16))ggsave(filename = "20240612_All_Prox_stacked_vln.pdf", plot = stack_vln, width = 18, height = 12, dpi = 600)#### Figure Supp 3d: Proximal Epi Cell Types ####epi_umap <- DimPlot(object = Epi_Named,                # Seurat object  reduction = 'umap',                 # Axes for the plot (UMAP, PCA, etc.) #group.by = "Patient",       # Labels to color the cells by ("seurat_clusters", "Age", "Time.Point)  repel = TRUE,                       # Whether to repel the cluster labelslabel = FALSE,                       # Whether to have cluster labels cols = c( "#35EFEF","#E95FE0","#B20224", "#F28D86", "#FB1111", "#FEB0DB"), pt.size = 1.6,                      # Size of each dot is (0.1 is the smallest)label.size = 2) +                   # Font size for labels    # You can add any ggplot2 1customizations herelabs(title = 'Colored by Cluster')+        # Plot titleNoLegend()ggsave(filename = "FIGs3b_epi_proximal_umap.pdf", plot = epi_umap, width = 15, height = 12, dpi = 600)#### Figure Supp 4e: Epi Markers All Distal Cells ####P_Epi_Filter <- readRDS(file = "../Proximal/20220914_Proximal_Epi_Cells.rds" , refhook =  NULL)P_Epi_Named <- RenameIdents(P_Epi_Filter, '0' = "Dbi+/Spdefhigh Secretory", '1' = "Bmpr1b+ Progenitor", '2' = "Wfdc2+ Secretory",'3' = "Ciliated", '4' = "Sox17high Secretory", '5' = "Kcne3+ Secretory")P_Epi_Named@active.ident <- factor(x = P_Epi_Named@active.ident, levels = c("Bmpr1b+ Progenitor","Wfdc2+ Secretory","Sox17high Secretory","Kcne3+ Secretory","Dbi+/Spdefhigh Secretory","Ciliated"))stack_vln <- StackedVlnPlot(obj = P_Epi_Named, features = features, slot = "data",pt.size = 0,cols = c( "#35EFEF","#F28D86", "#FB1111","#FEB0DB","#B20224","#E95FE0"))+theme(plot.title = element_text(size = 32, face = "bold.italic"))+scale_x_discrete(limits = c("Bmpr1b+ Progenitor","Wfdc2+ Secretory","Sox17high Secretory","Kcne3+ Secretory","Dbi+/Spdefhigh Secretory","Ciliated"))+theme(axis.text.x = element_text(size = 16, angle = 60))+theme(axis.text.y = element_text(size = 14))+theme(axis.title.y.left = element_text(size = 16))ggsave(filename = "20240612_Prox_Epi_stacked_vln.pdf", plot = stack_vln, width = 18, height = 12, dpi = 600)#### Figure Supp 4f: Proximal Epi Features####named_features <- c("Krt8","Epcam", "Msln","Slc1a3","Sox9","Itga6", "Bmpr1b","Ovgp1","Sox17","Pax8", "Egr1","Wfdc2","Dbi","Gsto1","Fxyd4","Vim","Kcne3","Spdef","Lgals1","Upk1a", "Thrsp","Selenop", "Gstm2","Anpep", "Klf6","Id2","Ifit1","Prom1", "Ly6a", "Kctd8", "Adam8","Foxj1","Fam183b","Rgs22","Dnali1", "Mt1" , "Dynlrb2")prox_dp <- DotPlot(object = Epi_Named,                    # Seurat objectassay = 'RNA',                        # Name of assay to use.  Default is the active assayfeatures = named_features,                  # List of features (select one from above or create a new one)# Colors to be used in the gradientcol.min = 0,                       # Minimum scaled average expression threshold (everything smaller will be set to this)col.max = 2.5,                        # Maximum scaled average expression threshold (everything larger will be set to this)dot.min = 0,                          # The fraction of cells at which to draw the smallest dot (default is 0)dot.scale = 6,                        # Scale the size of the pointsgroup.by = NULL,              # How the cells are going to be groupedsplit.by = NULL,                      # Whether to split the data (if you fo this make sure you have a different color for each variable)scale = TRUE,                         # Whether the data is scaledscale.by = "radius",                  # Scale the size of the points by 'size' or 'radius'scale.min = NA,                       # Set lower limit for scalingscale.max = NA                        # Set upper limit for scaling
)+    labs(x = NULL,                              # x-axis labely = NULL)+scale_color_viridis_c(option="F",begin=.4,end=0.9, direction = -1)+geom_point(aes(size=pct.exp), shape = 21, colour="black", stroke=0.6)+theme_linedraw()+guides(x =  guide_axis(angle = 90))+ theme(axis.text.x = element_text(size = 12 , face = "italic"))+theme(axis.text.y = element_text(size = 12))+theme(legend.title = element_text(size = 12))+ scale_y_discrete(limits = c("Ciliated","Dbi+/Spdefhigh Secretory","Kcne3+ Secretory","Sox17high Secretory","Wfdc2+ Secretory","Bmpr1b+ Progenitor"))ggsave(filename = "FIGs3c_epi_proximal_dotplot.pdf", plot = prox_dp, width = 18, height = 10, dpi = 600)

附图5

#### Figure Supp 5: Distal and proximal epithelial cell correlation ######## Packages Load ####
library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)#### Proximal Datasets ####Proximal <- readRDS( file = "../dataset/Proximal_Filtered_Cells.rds" , refhook =  NULL)Proximal_Named <- RenameIdents(Proximal, '0' = "Fibroblast 1", '1' = "Stem-like Epithelial", '2' = "Fibroblast 2",'3' = "Fibroblast 3", '4' = "Immune", '5' = "Secretory Epithelial", '6' = "Endothelial",'7' = "Ciliated Epithelial",'8' = "Mesothelial", '9' = "Smooth Muscle")Proximal_Named@active.ident <- factor(x = Proximal_Named@active.ident, levels = c('Fibroblast 1','Fibroblast 2','Fibroblast 3','Smooth Muscle','Endothelial','Stem-like Epithelial','Secretory Epithelial','Ciliated Epithelial','Immune','Mesothelial'))Proximal_Named <- SetIdent(Proximal_Named, value = Proximal_Named@active.ident)Epi_Filter <- readRDS(file = "../dataset/Proximal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Dbi+/Spdefhigh Secretory", '1' = "Bmpr1b+ Progenitor", '2' = "Wfdc2+ Secretory",'3' = "Ciliated", '4' = "Sox17high Secretory", '5' = "Kcne3+ Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c("Ciliated","Dbi+/Spdefhigh Secretory","Kcne3+ Secretory","Sox17high Secretory","Wfdc2+ Secretory","Bmpr1b+ Progenitor"))#### Proximal vs Distal Cluster Correlation ####Distal_Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Distal_Epi_Named <- RenameIdents(Distal_Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Distal_Epi_Named@active.ident <- factor(x = Distal_Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))prox_avg_exp <- AverageExpression(Epi_Named)$RNA
distal_avg_exp <- AverageExpression(Distal_Epi_Named)$RNAcor.exp <- as.data.frame(cor(x = prox_avg_exp , y = distal_avg_exp))cor.exp$x <- rownames(cor.exp)cor.df <- tidyr::gather(data = cor.exp, y, correlation, c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2"))distal_cells <- c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")prox_cells <- c("Bmpr1b+ Progenitor","Ciliated","Dbi+/Spdefhigh Secretory","Wfdc2+ Secretory","Sox17high Secretory","Kcne3+ Secretory")corr_matrix <- ggplot(cor.df, aes(x, y, fill = correlation)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_viridis_c(values = c(0,1),option="rocket", begin=.4,end=0.99, direction = -1,)+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 60, hjust = 1, size = 12, face = "bold"),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 12, face = "bold.italic"))+theme(plot.title = element_blank())+scale_y_discrete(limits = c("Ciliated 2","Ciliated 1","Selenop+/Gstm2high Secretory","Spdef+ Secretory","Fibroblast-like","Pax8low/Prom1+ Cilia-forming", "Slc1a3med/Sox9+ Cilia-forming","Cebpdhigh/Foxj1- Progenitor","Slc1a3+ Stem/Progenitor"))+scale_x_discrete(limits = c("Bmpr1b+ Progenitor","Wfdc2+ Secretory","Sox17high Secretory","Kcne3+ Secretory","Dbi+/Spdefhigh Secretory","Ciliated"))+geom_text(aes(x, y, label = round(correlation, digits = 2)), color = "black", size = 4)ggsave(filename = "FIGs3d_epi_cluster_corr.pdf", plot = corr_matrix, width = 18, height = 10, dpi = 600)

附图6

#### Figure Supp 6 and 10: Stem and Cancer Markers ######## Packages Load ####library(dplyr)
library(patchwork)
library(Seurat)
library(harmony)
library(ggplot2)
library(cowplot)
library(SoupX)
library(DoubletFinder)
library(data.table)
library(parallel)
library(tidyverse)
library(SoupX)
library(ggrepel)library(ggplot2)
library(gplots)
library(RColorBrewer)
library(viridisLite)
library(Polychrome)
library(circlize)
library(NatParksPalettes)library(monocle3)
library(ComplexHeatmap)
library(ggExtra)
library(gridExtra)
library(egg)library(scales)#### Distal Epithelial and Pseudotime Dataset ####Epi_Filter <- readRDS(file = "../dataset/Distal_Epi_Cells.rds" , refhook =  NULL)Epi_Named <- RenameIdents(Epi_Filter, '0' = "Spdef+ Secretory", '1' = "Slc1a3+ Stem/Progenitor", '2' = "Cebpdhigh/Foxj1- Progenitor",'3' = "Ciliated 1", '4' = "Ciliated 2", '5' = "Pax8low/Prom1+ Cilia-forming", '6' = "Fibroblast-like",'7' = "Slc1a3med/Sox9+ Cilia-forming",'8' = "Selenop+/Gstm2high Secretory")Epi_Named@active.ident <- factor(x = Epi_Named@active.ident, levels = c( c("Slc1a3+ Stem/Progenitor","Cebpdhigh/Foxj1- Progenitor","Slc1a3med/Sox9+ Cilia-forming","Pax8low/Prom1+ Cilia-forming", "Fibroblast-like","Spdef+ Secretory","Selenop+/Gstm2high Secretory","Ciliated 1","Ciliated 2")))cds <- readRDS(file = "../dataset/Distal_Epi_PHATE_Monocle3.rds" , refhook = NULL)#### Figure Supp 4: Stem Dot Plot ####stem_features <- c("Krt5","Krt17","Cd44","Prom1","Kit","Aldh1a1","Aldh1a2","Aldh1a3","Efnb1","Ephb1","Trp63","Sox2","Sox9","Klf4","Rnf43","Foxm1","Pax8","Nanog","Itga6","Psca","Tcf3","Tcf4","Nrp1","Slc1a3","Tnfrsf19","Smo","Lrig1","Ezh2","Egr1","Tacstd2","Dusp1","Slc38a2","Malat1","Btg2","Cdkn1c","Pdk4","Nedd9","Fos","Jun","Junb","Zfp36","Neat1","Gadd45g","Gadd45b")stem_dp <- DotPlot(object = Epi_Named,                    # Seurat objectassay = 'RNA',                        # Name of assay to use.  Default is the active assayfeatures = stem_features,                  # List of features (select one from above or create a new one)# Colors to be used in the gradientcol.min = 0,                       # Minimum scaled average expression threshold (everything smaller will be set to this)col.max = 2.5,                        # Maximum scaled average expression threshold (everything larger will be set to this)dot.min = 0,                          # The fraction of cells at which to draw the smallest dot (default is 0)dot.scale = 6,                        # Scale the size of the pointsgroup.by = NULL,              # How the cells are going to be groupedsplit.by = NULL,                      # Whether to split the data (if you fo this make sure you have a different color for each variable)scale = TRUE,                         # Whether the data is scaledscale.by = "radius",                  # Scale the size of the points by 'size' or 'radius'scale.min = NA,                       # Set lower limit for scalingscale.max = NA )+                       # Set upper limit for scalinglabs(x = NULL,                              # x-axis labely = NULL)+scale_color_viridis_c(option="F",begin=.4,end=0.9, direction = -1)+geom_point(aes(size=pct.exp), shape = 21, colour="black", stroke=0.6)+#theme_linedraw()+guides(x =  guide_axis(angle = 90))+theme(axis.text.x = element_text(size = 8 , face = "italic"))+theme(axis.text.y = element_text(size = 9))+theme(legend.title = element_text(size = 9))+theme(legend.text = element_text(size = 8))+ scale_y_discrete(limits = c("Ciliated 2","Ciliated 1","Selenop+/Gstm2high Secretory","Spdef+ Secretory","Fibroblast-like","Pax8low/Prom1+ Cilia-forming", "Slc1a3med/Sox9+ Cilia-forming","Cebpdhigh/Foxj1- Progenitor","Slc1a3+ Stem/Progenitor"))ggsave(filename = "FIGs4_stem_dp.pdf", plot = stem_dp, width = 12, height = 6, dpi = 600)x <- stem_dp$datawrite.csv( x , 'stem_dp_data.csv')#### Figure Supp 6: HGSC Driver Gene by Pseudotime ###### Calculate Pseudotime Values ##pseudo <- pseudotime(cds)Distal_PHATE@meta.data$Pseudotime <- pseudo # Add to Seurat Metadata## Subset Seurat Object ##color_cells <- DimPlot(Distal_PHATE , reduction = "phate", cols = c("#B20224", #1"#35EFEF", #2"#00A1C6", #3"#A374B5", #4"#9000C6", #5"#EA68E1", #6"lightgrey", #7"#2188F7", #8"#F28D86"),pt.size = 0.7,shuffle = TRUE,seed = 0,label = FALSE)## Psuedotime and Lineage Assignment ##cellID <- rownames(Distal_PHATE@reductions$phate@cell.embeddings)
phate_embeddings <- Distal_PHATE@reductions$phate@cell.embeddings
pseudotime_vals <- Distal_PHATE@meta.data$Pseudotimecombined_data <- data.frame(cellID, phate_embeddings, pseudotime_vals)# Calculate the Average PHATE_1 Value for Pseudotime Points = 0 #
avg_phate_1 <- mean(phate_embeddings[pseudotime_vals == 0, 1])# Pseudotime Values lower than avge PHATE_1 Embedding will be Negative to split lineages
combined_data$Split_Pseudo <- ifelse(phate_embeddings[, 1] < avg_phate_1, -pseudotime_vals, pseudotime_vals)# Define Lineage #
combined_data$lineage <- ifelse(combined_data$PHATE_1 < avg_phate_1, "Secretory",ifelse(combined_data$PHATE_1 > avg_phate_1, "Ciliogenic", "Progenitor"))Distal_PHATE$Pseudotime_Adj <- combined_data$Split_Pseudo
Distal_PHATE$Lineage <- combined_data$lineage# Subset #Pseudotime_Lineage <- subset(Distal_PHATE, idents = c("Secretory 1","Secretory 2","Msln+ Progenitor","Slc1a3+/Sox9+ Cilia-forming","Pax8+/Prom1+ Cilia-forming","Progenitor","Ciliated 1","Ciliated 2"))## Set Bins ##bins <- cut_number(Pseudotime_Lineage@meta.data$Pseudotime_Adj , 40) # Evenly distribute bins Pseudotime_Lineage@meta.data$Bin <- bins # Metadata for Bins## Set Idents to PSeudoime Bin ##time_ident <- SetIdent(Pseudotime_Lineage, value = Pseudotime_Lineage@meta.data$Bin)av.exp <- AverageExpression(time_ident, return.seurat = T)$RNA # Calculate Avg log normalized expression# Calculates Average Expression for Each Bin #
# if you set return.seurat=T, NormalizeData is called which by default performs log-normalization #
# Reported as avg log normalized expression ### Pseudotime Scale Bar ##list <- 1:40
colors = c(rev(rainbow20),rainbow20)
df <- data.frame(data = list, color = colors)pseudo_bar <- ggplot(df, aes(x = 1:40, y = 1, fill = color)) + geom_bar(stat = "identity",position = "fill", color = "black", size = 0, width = 1) +scale_fill_identity() +theme_void()+ theme(axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())ggsave(filename = "pseudo_bar.pdf", plot = pseudo_bar, width = 0.98, height = 0.19, dpi = 600)## Plot HGSC driver gene list across pseudotime bin ##features <- c("Trp53", "Brca1", "Brca2",	"Csmd3",	"Nf1",	"Fat3",	"Gabra6", "Rb1", "Apc",	"Lrp1b","Prim2",	"Cdkn2a", "Crebbp",	"Wwox", "Ankrd11",	"Map2k4",	"Fancm",	"Fancd2",	"Rad51c",  "Pten")# Create Bin List and expression of features #bin_list <- unique(Pseudotime_Lineage@meta.data$Bin) plot_info <- as.data.frame(av.exp[features,]) # Call Avg Expression for featuresz_score <- transform(plot_info, SD=apply(plot_info,1, mean, na.rm = TRUE))
z_score <- transform(z_score, MEAN=apply(plot_info,1, sd, na.rm = TRUE))z_score1 <- (plot_info-z_score$MEAN)/z_score$SDplot_info$y <- rownames(plot_info) # y values as features
z_score1$y <- rownames(plot_info)plot_info <- gather(data = plot_info, x, expression, bin_list) #set plot
z_score1 <- gather(data = z_score1, x, z_score, bin_list) #set plot# Create Cell Clusters DF #Labeled_Pseudotime_Lineage <- RenameIdents(Pseudotime_Lineage, 'Secretory 1' = "Spdef+ Secretory", 'Progenitor' = "Slc1a3+ Stem/Progenitor", 'Msln+ Progenitor' = "Cebpdhigh/Foxj1- Progenitor",'Ciliated 1' = "Ciliated 1", 'Ciliated 2' = "Ciliated 2", 'Pax8+/Prom1+ Cilia-forming' = "Pax8low/Prom1+ Cilia-forming", 'Fibroblast-like' = "Fibroblast-like", #removed'Slc1a3+/Sox9+ Cilia-forming' = "Slc1a3med/Sox9+ Cilia-forming",'Secretory 2' = "Selenop+/Gstm2high Secretory")cluster_table <- table(Labeled_Pseudotime_Lineage@active.ident, Labeled_Pseudotime_Lineage@meta.data$Bin)clusters <- data.frame(cluster_table)clusters <- clusters %>% group_by(Var2) %>%mutate(Perc = Freq / sum(Freq))# Create Pseudotime DF #pseudotime_table <- table(seq(1, length(bin_list), 1), unique(Labeled_Pseudotime_Lineage@meta.data$Bin),seq(1, length(bin_list), 1))pseudotime_bins <- data.frame(pseudotime_table)  # calculate max and min z-scores
max_z <- max(z_score1$z_score, na.rm = TRUE)
min_z <- min(z_score1$z_score, na.rm = TRUE)# set color for outliers
outlier_color <- ifelse(z_score1$z_score > max_z | z_score1$z_score < min_z, ifelse(z_score1$z_score > 0, "#AD1F24", "#51A6DC"), "#e2e2e2")## Plot Gene Expression ### Set different na.value options for positive and negative values
na_color_pos <- "#AD1F24" # color for positive NA values
na_color_neg <- "#51A6DC" # color for negative NA valuescustom_bin_names <- c(paste0("S", 20:1), paste0("C", 1:20))figure <- ggplot(z_score1, aes(x, y, fill = z_score)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradientn(colors=c("#1984c5", "#e2e2e2", "#c23728"), name = "Average Expression \nZ-Score", limits = c(-3,3), na.value = ifelse(is.na(z_score1) & z_score1 > 0, na_color_pos, ifelse(is.na(z_score1) & z_score1 < 0, na_color_neg, "grey50")),oob = scales::squish)+scale_x_discrete(limits= sort(bin_list) , labels= custom_bin_names)+scale_y_discrete(limits= rev(features))+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),                                      # Text size throughout the plotaxis.text.x = element_text(color = 'black', angle = 0, hjust = 0.5, size = 10, face = "bold"),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold.italic"))+theme(plot.title = element_blank(),plot.margin=unit(c(-0.5,1,1,1), "cm"))## Plot Cluster Percentage ##`Spdef+ Secretory` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Spdef+ Secretory")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(1,1,1,1), "cm"))`Selenop+/Gstm2high Secretory` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Selenop+/Gstm2high Secretory")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Cebpdhigh/Foxj1- Progenitor` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Cebpdhigh/Foxj1- Progenitor")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Slc1a3+ Stem/Progenitor` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Slc1a3+ Stem/Progenitor")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Slc1a3med/Sox9+ Cilia-forming` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Slc1a3med/Sox9+ Cilia-forming")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Pax8low/Prom1+ Cilia-forming` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Pax8low/Prom1+ Cilia-forming")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Ciliated 1` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Ciliated 1")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))`Ciliated 2` <- ggplot(clusters, aes(Var2, Var1, fill = Perc)) +geom_tile(color = "black",lwd = 1,linetype = 1) +scale_fill_gradient2(low="white", high="#000000", mid = "white", midpoint = 0, name = "Percentage")+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Ciliated 2")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust = 1, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))## Plot Pseudotime Color ##list <- 1:40
colors = c(rev(rainbow20),rainbow20)
df <- data.frame(data = list, color = colors)binning <- ggplot(df, aes(x = 1:40, y = 1, fill = color)) + geom_bar(stat = "identity",position = "fill", color = "black", size = 1, width = 1) +scale_fill_identity() +theme_void()+ theme(axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())+scale_x_discrete(limits= sort(bin_list) , labels= seq(1, length(bin_list), 1))+scale_y_discrete(limits= "Pseudotime Bin ")+theme(panel.background = element_blank())+labs(title = "Expression of Genes by Pseudotime Bin" ,x = element_blank(),y = element_blank())+theme(text = element_text(size = 12, face = "bold"),# Text size throughout the plotaxis.ticks.x = element_blank(),axis.ticks.y = element_blank(),axis.text.x = element_blank(),    # Text color, angle, and horizontal adjustment on x-axis axis.text.y = element_text(color = 'black', hjust =1, vjust = .75, size = 14, face = "bold"))+theme(plot.title = element_blank(),plot.margin=unit(c(-1.25,1,1,1), "cm"))### Combine Plots ###psuedotime_lineage <- ggarrange(`Spdef+ Secretory`,`Selenop+/Gstm2high Secretory`,`Cebpdhigh/Foxj1- Progenitor`,`Slc1a3+ Stem/Progenitor`,`Slc1a3med/Sox9+ Cilia-forming`,`Pax8low/Prom1+ Cilia-forming`,`Ciliated 1`,`Ciliated 2`,`binning`,figure , ncol=1,heights = c(2, 2, 2, 2, 2, 2, 2, 2, 2, (2*length(features)),widths = c(3)),padding = unit(0.01))ggsave(filename = "FIGs6_psuedotime_driver_gene.pdf", plot = psuedotime_lineage, width = 18, height = 9, dpi = 600)write.csv(z_score1 , 'cancer_pseudotime.csv')