-
Notifications
You must be signed in to change notification settings - Fork 6
/
Figure E07
150 lines (136 loc) · 15.7 KB
/
Figure E07
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#Object
##############################################################################################################################################################################
library(monocle3); library(monocle); library(SeuratWrappers); library(CytoTRACE);
Original.OView03<- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_OverView03.RDS");
CD4THluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_CD4THluster01.RDS");
CD8THluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_CD8THluster01.RDS");
NHluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_NKHluster03.RDS");
BHluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_IntBHluster04_RMIG.RDS");
MHluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_MyeloidHluster01.RDS");
StrHluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_StrHluster01.RDS");
##############################################################################################################################################################################
#Function
##############################################################################################################################################################################
PaperBox = function(DataFrame, Color, Method, Title)
{
Plot <- ggplot(DataFrame, aes(Class, Value)) + geom_boxplot(aes(color=Class), notch=FALSE, outlier.shape=NA, width=0.5) + geom_jitter(aes(colour = Class), size = 1.5, alpha=0.8, width=0.3) +
scale_fill_manual(values=Color) + theme_classic() + stat_compare_means(method = Method) + ggtitle(Title) + # + labs(x=Axis[1], y=Axis[2])
theme(axis.text = element_text(size=10,face="plain"), axis.text.x = element_text(angle=45, hjust=1), axis.text.y = element_blank(), #axis.title=element_text(size=14,face="plain",color="black"),
plot.title = element_text(hjust = 0.5, size=8, face="plain"), axis.title = element_blank(), legend.position="none")
return(Plot);
}
MultiBox = function(DataFrame, GroupColor, NodeColor, LineColor, Title, Method)
{
Plot <- ggplot(DataFrame, aes(TissueType, Value)) + geom_boxplot(aes(color=TissueType), color=GroupColor, notch=FALSE, outlier.shape=NA, width=0.5) + geom_point(aes(colour = TissueLabl), color=NodeColor, size = 1.5, alpha=0.8, width=0.3) +
theme_classic() + stat_compare_means(method = Method) + geom_line(aes(group=Patient, colour="#d9d9d9"), color="#d9d9d9", linetype="11") + ggtitle(Title) + # + labs(x=Axis[1], y=Axis[2])
theme(axis.text = element_text(size=10,face="plain"), axis.text.y = element_blank(), #axis.title=element_text(size=14,face="plain",color="black"), axis.text.x = element_text(angle=45, hjust=1),
plot.title = element_text(hjust = 0.5, size=8, face="plain"), axis.title = element_blank(), legend.position="none")
return(Plot);
}
##############################################################################################################################################################################
#Fig E7a
##############################################################################################################################################################################
NormalHluster <- AddInfo(NormalHluster); CancerCluster <- AddInfo(CancerCluster); #DimPlot for Non-Malignan Cells
Color <- c("#fdbf6f","#cab2d6","#6a3d9a","#1f78b4","#ffff99","#e31a1c","#b2df8a","#a6cee3","#fb9a99","#b15928","#33a02c","#ff7f00");
Cluster <- c("Quiet_Foveolar","Parietal","MN/BG","Chief","Endocrine","Tuft","Quiet_KRT7+","Goblet","Enterocyte","Hepatocyte","Prolif_KRT7+","Prolif_Foveolar");
HPpositive <- rownames(NormalHluster[[]])[which(NormalHluster[[]]$HPinfec=="Positive")]; HPnegative <- rownames(NormalHluster[[]])[which(NormalHluster[[]]$HPinfec=="Negative")];
pdf("/rsrch3/scratch/genomic_med/edai/Project004/Figure/Lineage_01.pdf", width=10, height=10)
Idents(NormalHluster) <- "CellType"; Idents(NormalHluster) <- factor(Idents(NormalHluster), levels=Cluster);
Plot01 <- DimPlot(object=NormalHluster, reduction="umap", cols=Color, pt.size=1.0, label=F) + xlim(-8,8) + theme_void() + #ylim(-15,17.5) +
theme(legend.position="none",axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(),
axis.text.y = element_blank());
Combine <- plot_grid(Plot01,nrow=1,ncol=1); print(Combine);
dev.off()
##############################################################################################################################################################################
#Fig E7b
##############################################################################################################################################################################
Idents(NormalHluster) <- "CellType"; Objt <- subset(NormalHluster,idents=c("Quiet_KRT7+","Prolif_KRT7+","Prolif_Foveolar","Quiet_Foveolar"));
saveRDS(Objt,"/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_TempObjt.RDS")
Objt <- NormalizeData(Objt, normalization.method = "LogNormalize", scale.factor = 10000); Objt <- ScaleData(Objt, features=rownames(Objt), verbose = FALSE);
Idents(Objt) <- "CellType"; Idents(Objt) <- factor(Idents(Objt), levels=c("Quiet_KRT7+","Prolif_KRT7+","Prolif_Foveolar","Quiet_Foveolar"));
EpiMarkers <- FindAllMarkers(Objt, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25); EpiMarkers %>% group_by(cluster) %>% top_n(n = 30, wt = avg_log2FC) -> Top30;
#write.table(EpiMarkers, file="/rsrch3/scratch/genomic_med/edai/Project004/Temporary/MarkerList.csv", sep=",", row.names=T, col.names=T);
pdf("/rsrch3/scratch/genomic_med/edai/Project004/Figure/Feature_01.pdf", width=25, height=15)
DoHeatmap(Objt,label = TRUE, features = Top30$gene) + NoLegend()
dev.off()
##############################################################################################################################################################################
#Trajectory Analysis
##############################################################################################################################################################################
Objt <- list(); Idents(NormalHluster) <- "CellType"; Objt[[1]] <- subset(NormalHluster, idents=c("Quiet_Foveolar","Prolif_Foveolar","Quiet_KRT7+","Prolif_KRT7+"));
Idents(Objt[[1]]) <- "TissueType"; Objt[[2]] <- subset(Objt[[1]], idents=c("Adjacent")); Objt[[3]] <- subset(Objt[[1]], idents=c("Primary")); Objt[[4]] <- subset(Objt[[1]], idents=c("Metastasis"));
Objt[[5]] <- subset(Objt[[1]], idents=c("Adjacent","Primary")); Objt[[6]] <- subset(Objt[[1]], idents=c("Primary","Metastasis")); names(Objt) <- c("Ad&Pr&Mt","Ad","Pr","Mt","Ad&Pr","Pr&Mt");
saveRDS(ObjtList, "/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_Fig1Trajectory.RDS");
ObjtList <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_Fig1Trajectory.RDS"); Objt <- ObjtList[[2]]; ExprMatrix <- Objt@assays$RNA@counts; SampleInfo <- Objt@meta.data;
GeneAnno <- data.frame(gene_short_name=rownames(Objt)); rownames(GeneAnno) <- GeneAnno$gene_short_name; PD <- new("AnnotatedDataFrame", data = SampleInfo); FD <- new("AnnotatedDataFrame", data = GeneAnno);
CellDataSetV2 <- newCellDataSet(ExprMatrix, phenoData = PD, featureData = FD, expressionFamily = negbinomial.size()); CellDataSetV2 <- estimateSizeFactors(CellDataSetV2); CellDataSetV2 <- estimateDispersions(CellDataSetV2);
CellDataSetV2 <- setOrderingFilter(CellDataSetV2, ordering_genes = VariableFeatures(Objt)); CellDataSetV2 <- reduceDimension(CellDataSetV2, method = 'DDRTree'); CellDataSetV2 <- orderCells(CellDataSetV2);
CellDataSetV3 <- as.cell_data_set(Objt); CellDataSetV3 <- cluster_cells(CellDataSetV3, reduction_method = "UMAP"); CellDataSetV3 <- learn_graph(CellDataSetV3, use_partition = TRUE);
Idents(Objt) <- "CellType"; Root <- WhichCells(Objt, idents="Quiet_KRT7+"); CellDataSetV3 <- order_cells(CellDataSetV3, reduction_method = "UMAP", root_cells = Root);
#diff_test_res <- differentialGeneTest(CellDataSetV2); save(Objt,CellDataSetV2,CellDataSetV3,diff_test_res,file="/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Figure1_Trajectory.RDA")
GeneExpression <- as.matrix(GetAssayData(Objt)); CytoTRACE <- CytoTRACE(GeneExpression);
Mono3.Pseudotime <- pseudotime(CellDataSetV3, reduction_method = "UMAP"); phenoData(CellDataSetV2)$Monocle3 <- Mono3.Pseudotime; phenoData(CellDataSetV2)$CytoTRACE <- CytoTRACE[[1]];
Trajectory <- data.frame(Monocle2=phenoData(CellDataSetV2)$Pseudotime,Monocle3=phenoData(CellDataSetV2)$Monocle3,CytoTRACE=phenoData(CellDataSetV2)$CytoTRACE); rownames(Trajectory) <- rownames(Objt[[]])
save(Objt,CellDataSetV2,CellDataSetV3,CytoTRACE,Trajectory,file="/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_Fig1Trajectory.RDA");
##############################################################################################################################################################################
#Fig E7c
##############################################################################################################################################################################
#left
load("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_Fig1Trajectory.RDA"); Objt <- AddMetaData(Objt,Trajectory,col.name=colnames(Trajectory));
pdf("/rsrch3/scratch/genomic_med/edai/Project004/Figure/Feature_01.pdf", width=12, height=4)
Plot01 <- plot_cell_trajectory(CellDataSetV2, color_by = "CellType", cell_size=0.1) + geom_point(aes(color = CellType), size=0.0001) +
scale_color_manual(breaks = c("Quiet_Foveolar","Quiet_KRT7+","Prolif_KRT7+","Prolif_Foveolar"),
values=c("#fdbf6f","#b2df8a","#33a02c","#ff7f00")) + theme(legend.position="right");
Plot02 <- plot_cell_trajectory(CellDataSetV2, color_by = "Monocle2", cell_size=0.1) + scale_colour_gradientn(colours = rev(rainbow(100)[c(1:70)])) + theme(legend.position="right");
Plot03 <- plot_cell_trajectory(CellDataSetV2, color_by = "Monocle3", cell_size=0.1) + scale_colour_gradientn(colours = rev(rainbow(100)[c(1:70)])) + theme(legend.position="right");
Plot04 <- plot_cell_trajectory(CellDataSetV2, color_by = "CytoTRACE", cell_size=0.1) + scale_colour_gradientn(colours = rainbow(100)[c(1:70)]) + theme(legend.position="right");
Combine <- plot_grid(Plot01,Plot02,Plot03,Plot04,nrow=1,ncol=4); print(Combine);
dev.off()
#right
SmoothPlot = function(DataFrame,Title)
{
Plot <- ggplot(data = DataFrame, aes(X,Y)) + geom_point(fill="black",colour="black",size=0.5,shape=21) + geom_smooth(method = 'loess',span=0.1,se=TRUE,colour="#00A5FF",fill="#00A5FF",alpha=0.2) +
scale_y_continuous(breaks = seq(0, 125, 25)) + ggtitle(Title) + theme_classic() +
theme(panel.background=element_rect(fill="white",colour="white",size=0), axis.line=element_line(colour="black",size=0.25), axis.title=element_text(size=14,face="plain",color="black"),
axis.text = element_text(size=10,face="plain"), axis.title.y = element_blank(), axis.title.x = element_blank(), legend.position="none");
return(Plot)
}
SmoothMerge = function(DataFrame,Title,Color)
{
Plot <- ggplot(DataFrame, aes(x,value,group=variable,col=variable)) + geom_smooth(method = 'loess', span=0.5, se=TRUE, size=1) +
scale_colour_manual(values=Color) + ggtitle(Title) +
theme(panel.background=element_rect(fill="white",colour="white",size=0), axis.line=element_line(colour="black",size=0.25),
axis.title=element_blank(),axis.text = element_blank(), legend.position="right"); return(Plot);
};
load("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_Fig1Trajectory.RDA"); #Objt <- ScaleData(Objt, features=SelectMark, verbose = FALSE);
SelectMark <- c("OLFM4","ALDH1A1","TFF2","MUC5AC","GKN1","GKN2"); ScaleMatrix <- t(as.matrix(GetAssayData(Objt, slot = "data")));
ScaleMatrix <- ScaleMatrix[match(rownames(Trajectory),rownames(ScaleMatrix)),match(SelectMark,colnames(ScaleMatrix))]; Trajectory <- cbind(Trajectory,as.data.frame(ScaleMatrix))
Color <- c("OLFM4"="#66c2a5","ALDH1A1"="#fc8d62","TFF2"="#8da0cb","MUC5AC"="#e78ac3")
PlotData01 <- data.frame(x=Trajectory$Monocle2,value=Trajectory$OLFM4,variable=rep("OLFM4",nrow(Trajectory)));
PlotData02 <- data.frame(x=Trajectory$Monocle2,value=Trajectory$ALDH1A1,variable=rep("ALDH1A1",nrow(Trajectory)));
PlotData03 <- data.frame(x=Trajectory$Monocle2,value=Trajectory$TFF2,variable=rep("TFF2",nrow(Trajectory)));
PlotData04 <- data.frame(x=Trajectory$Monocle2,value=Trajectory$MUC5AC,variable=rep("MUC5AC",nrow(Trajectory)));
PlotData05 <- data.frame(x=Trajectory$CytoTRACE,value=Trajectory$OLFM4,variable=rep("OLFM4",nrow(Trajectory)));
PlotData06 <- data.frame(x=Trajectory$CytoTRACE,value=Trajectory$ALDH1A1,variable=rep("ALDH1A1",nrow(Trajectory)));
PlotData07 <- data.frame(x=Trajectory$CytoTRACE,value=Trajectory$TFF2,variable=rep("TFF2",nrow(Trajectory)));
PlotData08 <- data.frame(x=Trajectory$CytoTRACE,value=Trajectory$MUC5AC,variable=rep("MUC5AC",nrow(Trajectory)));
PlotData1_1 <- rbind(PlotData01,PlotData02); PlotObjt1_1 <- SmoothMerge(PlotData1_1,"Monocle2",c("OLFM4"="#66c2a5","ALDH1A1"="#fc8d62"));
PlotData1_2 <- rbind(PlotData03,PlotData04); PlotObjt1_2 <- SmoothMerge(PlotData1_2,"Monocle2",c("TFF2"="#8da0cb","MUC5AC"="#e78ac3"));
PlotData2_1 <- rbind(PlotData05,PlotData06); PlotObjt2_1 <- SmoothMerge(PlotData2_1,"CytoTRACE",c("OLFM4"="#66c2a5","ALDH1A1"="#fc8d62"));
PlotData2_2 <- rbind(PlotData07,PlotData08); PlotObjt2_2 <- SmoothMerge(PlotData2_2,"CytoTRACE",c("TFF2"="#8da0cb","MUC5AC"="#e78ac3"));
pdf("/rsrch3/scratch/genomic_med/edai/Project004/Figure/Original_01.pdf", width=8, height=6)
Combine <- plot_grid(PlotObjt1_1,PlotObjt1_2, PlotObjt2_1,PlotObjt2_2, nrow=2,ncol=2); print(Combine);
dev.off()
##############################################################################################################################################################################
#Fig E7d
##############################################################################################################################################################################
NormalHluster <- readRDS("/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/Original_HlusterNormal05.RDS"); Idents(NormalHluster) <- "CellType";
Objt <- subset(NormalHluster, idents=c("Quiet_Foveolar","Prolif_Foveolar","Quiet_KRT7+","Prolif_KRT7+")); ExprMatrix <- Objt@assays$RNA@counts; SampleInfo <- Objt@meta.data;
GeneAnno <- data.frame(gene_short_name=rownames(Objt)); rownames(GeneAnno) <- GeneAnno$gene_short_name; PD <- new("AnnotatedDataFrame", data = SampleInfo); FD <- new("AnnotatedDataFrame", data = GeneAnno);
CellDataSetV2 <- newCellDataSet(ExprMatrix, phenoData = PD, featureData = FD, expressionFamily = negbinomial.size()); CellDataSetV2 <- estimateSizeFactors(CellDataSetV2); CellDataSetV2 <- estimateDispersions(CellDataSetV2);
CellDataSetV2 <- setOrderingFilter(CellDataSetV2, ordering_genes = VariableFeatures(Objt)); CellDataSetV2 <- reduceDimension(CellDataSetV2, method = 'DDRTree'); CellDataSetV2 <- orderCells(CellDataSetV2);
diff_test_res <- differentialGeneTest(CellDataSetV2); save(CellDataSetV2,diff_test_res,file="/rsrch3/scratch/genomic_med/edai/Project004/DataSpace/FigureS4_Trajectory.RDA");
diff_test_res <- diff_test_res[which(diff_test_res$status=="OK"&diff_test_res$use_for_ordering==T&diff_test_res$qval<=1e-60),]; SigGene <- diff_test_res$gene_short_name;
pdf("/rsrch3/scratch/genomic_med/edai/Project004/Figure/Original_01.pdf", width=6, height=6)
plot_pseudotime_heatmap(CellDataSetV2[SigGene,], num_clusters = 4, cores = 1, show_rownames = F,return_heatmap = F);
dev.off()
##############################################################################################################################################################################