plot27 · 代码 20 · 图 0 · 笔记 3 · 代码库专题
CRC_scRNAseq
Single-cell transcriptomic profiling of human colorectal cancer, including the annotation of cell types, assessing the robustness of cell types classification, the function enrichment of cell types, i
CRC_scRNAseq（自动生成占位图）
图源：（自动生成占位图）
分析代码
8 个源码文件，按目录顺序；本机绝对路径已替换为 ~。 代码可复制，但需要按自己的数据改路径与参数。
1.maintype.r
WholeTissueDataN <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Process/WholeTissueData_Filter_nFeatureMito.rds.gz")
EpithelialCells <- c("EPCAM","KRT8","KRT18")
StromalCells <- c("COL1A1","COL1A2","COL6A1","COL6A2","VWF","PLVAP","CDH5","S100B")
ImmuneCells <- c("CD52","CD2","CD3D","CD3G","CD3E","CD79A","CD79B","CD14","FCGR3A","CD68","CD83","CSF1R","FCER1G")
MyeloidCells <- c("CD68", "XCR1", "CLEC9A", "CLEC10A", "CD1C", "S100A8", "S100A9", "TPSAB1", "OSM")
TCells <- c("NKG7", "KLRC1", "CCR7", "FOXP3", "CTLA4", "CD8B", "CXCR6", "CD3D")
BCells<- c("MZB1", "IGHA1", "SELL", "CD19", "AICDA")
Samples <- c("p7_N","p9_N", "p10_N", "p11_N", "p12_N", "p7_T", "p9_T", "p10_T", "p11_T", "p12_T")
WholeTissueDataN <- NormalizeData(WholeTissueDataN)
WholeTissueDataN <- FindVariableFeatures(object =WholeTissueDataN, mean.function = ExpMean, dispersion.function = LogVMR)
WholeTissueDataN <- ScaleData(object = WholeTissueDataN)
WholeTissueDataN <- RunPCA(object = WholeTissueDataN, npcs = 30, verbose = FALSE)
Idents(WholeTissueDataN) <- WholeTissueDataN$Pat_Tissues
ColonData <- WholeTissueDataN %>%
RunHarmony("Pat_Tissues", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(ColonData, 'harmony')
ColonData <- ColonData %>%
RunUMAP(reduction = "harmony", dims = 1:20) %>%
FindNeighbors(reduction = "harmony", dims = 1:20)
ColonData <- ColonData %>%
RunTSNE(reduction = "harmony", dims = 1:20)
for(i in seq(0.1,0.8,by=0.1)){
ColonData <- FindClusters(ColonData, resolution =i, algorithm = 1)
}
ColonData@meta.data$MyeloidCells<- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% MyeloidCells,],2,mean)
ColonData@meta.data$TCells <- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% TCells,],2,mean)
ColonData@meta.data$BCells <- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% BCells,],2,mean)
ColonData@meta.data$EpithelialCells <- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% EpithelialCells,],2,mean)
ColonData@meta.data$StromalCells<- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% StromalCells,],2,mean)
ColonData@meta.data$ImmuneCellsScore <- apply(ColonData@assays$RNA@data[rownames(ColonData@assays$RNA@data) %in% ImmuneCells,],2,mean)
ColonData@meta.data$ImmuneCells <- apply(ColonData@meta.data[,c("ImmuneCellsScore","MyeloidCells","TCells","BCells")],1,max)
readr::write_rds(ColonData,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColonData_Harmony.rds.gz",compress = "gz")
Types <- c("MyeloidCells","TCells", "BCells","ImmuneCellsScore","EpithelialCells", "StromalCells")
ThreeTypes <- c("ImmuneCells","EpithelialCells", "StromalCells")
ImmTypes <- c("MyeloidCells","TCells", "BCells")
p1 <- DimPlot( ColonData, reduction = "umap",label=T,group.by = "RNA_snn_res.0.8",
cols = colorRampPalette(.cluster_cols,space="rgb")(length(unique(ColonData$RNA_snn_res.0.8))))+NoLegend()
p <-FeaturePlot(ColonData,reduction = "umap",features = Types,cols = c("gray","orange","orange","red","red"),combine = FALSE)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend() + NoAxes()
}
# pdf(paste("Fig/",Sample,"Cluster_Score.pdf",sep=""),width = 8,height =8)
# cowplot::plot_grid(plotlist = p)
# dev.off()
# #cowplot::plot_grid(plotlist = p)
# pdf(paste("Fig/",Sample,"Cluster_VlnPlot.pdf",sep=""),width = 10,height =3)
ColonDataS$RNA_snn_res.0.8 <- factor(ColonDataS$RNA_snn_res.0.8,levels = 0:(length(unique(ColonDataS$RNA_snn_res.0.8))-1))
Vln_p <- VlnPlot(ColonData, features= Types,group.by = "RNA_snn_res.0.8",
pt.size = 0,cols=.cluster_cols[1:length(unique(ColonData@meta.data$RNA_snn_res.0.8))],combine = F)
# dev.off()
ClusterSum <- sapply(split(ColonData@meta.data[,ThreeTypes],ColonData@meta.data$RNA_snn_res.0.8),function(x){apply(x,2,mean)}) %>% t %>% data.frame %>%
dplyr::mutate(Cluster=rownames(.))
ClusterSum[,1:3] <- apply(ClusterSum[,1:3],2,function(x){(x)})
ClusterSum$ClusterBelong <- apply(ClusterSum[,1:3],1,function(x)match(max(x),x))
ClusterSum$ClusterDefine <- ifelse(ClusterSum$ClusterBelong==1,"ImmuneCells",ifelse(ClusterSum$ClusterBelong==2,"EpithelialCells","StromalCells"))
ClusterSum <- ClusterSum[order(ClusterSum$ClusterDefine),]
ClusterSum.m <- ClusterSum %>% tidyr::gather(key=ThreeTypes,value=ClusterScore,-c(Cluster,ClusterBelong,ClusterDefine))
ClusterDefineP <- ggplot(ClusterSum.m,aes(y=Cluster,x=ThreeTypes))+
geom_tile(aes(fill=ClusterScore))+
scale_x_discrete(limit=ThreeTypes)+
scale_y_discrete(limit=ClusterSum$Cluster)+
#coord_fixed(ratio = 2)+
scale_fill_gradientn(colors = colorRampPalette(c("gray","orange","red"),space = "rgb")(20),guide = F)+
theme(panel.background=element_blank(),
panel.grid=element_blank(),
axis.title=element_blank(),
axis.text.x=element_text(size=8,colour = "black",angle=90,hjust=1,vjust=0.5),
axis.text.y=element_text(size=8),
axis.ticks=element_blank())
##Immune types
ImmClusterSum <- sapply(split(ColonData@meta.data[,ImmTypes],ColonData@meta.data$RNA_snn_res.0.8),function(x){apply(x,2,mean)}) %>% t %>% data.frame %>%
dplyr::mutate(Cluster=rownames(.))
ImmClusterSum[,1:3] <- apply(ImmClusterSum[,1:3],2,function(x){(x)})
ImmClusterSum$ClusterBelong <- apply(ImmClusterSum[,1:3],1,function(x)match(max(x),x))
ImmClusterSum$ClusterDefine <- ClusterSum$ClusterDefine[match(ImmClusterSum$Cluster,ClusterSum$Cluster)]
ImmClusterSum$ClusterDefine <- ifelse(ImmClusterSum$ClusterDefine %in% c("EpithelialCells","StromalCells"),
ImmClusterSum$ClusterDefine,
ifelse(ImmClusterSum$ClusterBelong==1,"MyeloidCells",ifelse(ImmClusterSum$ClusterBelong==2,"TCells","BCells")))
ImmClusterSumM <- ImmClusterSum
ImmClusterSum <- ImmClusterSum[order(ImmClusterSum$ClusterDefine),]
ImmClusterSum <- ImmClusterSum[ImmClusterSum$ClusterDefine %in% ImmTypes,]
ImmClusterSum.m <- ImmClusterSum %>% tidyr::gather(key=ImmTypes,value=ClusterScore,-c(Cluster,ClusterBelong,ClusterDefine))
ImmTypesP <- ggplot(ImmClusterSum.m,aes(y=Cluster,x=ImmTypes))+
geom_tile(aes(fill=ClusterScore))+
scale_x_discrete(limit=ImmTypes)+
scale_y_discrete(limit=ImmClusterSum$Cluster)+
#coord_fixed(ratio = 2)+
scale_fill_gradientn(colors = colorRampPalette(c("gray","orange","red"),space = "rgb")(20),guide = F)+
theme(panel.background=element_blank(),
panel.grid=element_blank(),
axis.title=element_blank(),
axis.text.x=element_text(size=8,colour = "black",angle=90,hjust=1,vjust=0.5),
axis.text.y=element_text(size=8),
axis.ticks=element_blank())
Pre_l <- length(p)
for(i in 1:length(Vln_p)){
p[[(i+Pre_l)]] <- Vln_p[[i]] + NoLegend()+theme(axis.text.x = element_text(size=7))
}
p[[(length(p)+1)]] <- p1+ NoLegend() + NoAxes()
p[[(length(p)+1)]] <- ClusterDefineP
p[[(length(p)+1)]] <- ImmTypesP
###SingleR
ref.se <- readr::read_rds("~/Project/Collaboration/2020ShenLei/Matrix/ProcessData/SingeR_ref.se.rds.gz")
bp.se <- ref.se$BP
hpca.se <- ref.se$HPCA
ColonDataS <- ColonData
###PBMC
ColonDataS_SingleR_Main <- SingleR::SingleR(GetAssayData(ColonDataS@assays$RNA, slot = "data"),
ref=ref.se,
method = "cluster",
clusters =ColonDataS@meta.data$RNA_snn_res.0.8,
labels = list(bp.se$label.main, hpca.se$label.main))
SingleR_Result <- data.frame(HPCA =ColonDataS_SingleR_Main@listData$orig.results$HPCA$first.labels,
BP =ColonDataS_SingleR_Main@listData$orig.results$BP$first.labels,
clusters=sort(unique(ColonDataS@meta.data$RNA_snn_res.0.8)))
ColonDataS@meta.data$HPCA <- SingleR_Result$HPCA[match(ColonDataS@meta.data$RNA_snn_res.0.8,SingleR_Result$clusters)]
ColonDataS@meta.data$BP <- SingleR_Result$BP[match(ColonDataS@meta.data$RNA_snn_res.0.8,SingleR_Result$clusters)]
HPCA_p <- DimPlot( ColonDataS, reduction = "umap",label=T,group.by = "HPCA",cols = .cluster_cols)+NoLegend()
BP_p <- DimPlot( ColonDataS, reduction = "umap",label=T,group.by = "BP",cols = .cluster_cols)+NoLegend()
p[[(length(p)+1)]] <- HPCA_p
p[[(length(p)+1)]] <- BP_p
#####
ColonData@meta.data$MainTypes <- ClusterSum$ClusterDefine[match(ColonData@meta.data$RNA_snn_res.0.8,ClusterSum$Cluster)]
ColonData@meta.data$AllTypes <- ImmClusterSumM$ClusterDefine[match(ColonData@meta.data$RNA_snn_res.0.8,ImmClusterSumM$Cluster)]
MainType_p <- DimPlot( ColonData, reduction = "umap",label=T,group.by = "MainTypes",cols = .cluster_cols)+NoLegend()
AllType_p <- DimPlot( ColonData, reduction = "umap",label=T,group.by = "AllTypes",cols = .cluster_cols[c(9,1,12,3,8)])+NoLegend()
p[[(length(p)+1)]] <- MainType_p
p[[(length(p)+1)]] <- AllType_p
png("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/res0.8_Cluster.png",width = 4200,height =2100,res=200)
ggp <- cowplot::plot_grid(plotlist = p,ncol = 6)
print(ggp)
dev.off()
readr::write_rds(ColonData,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColonData_MaintypeDefine.rds.gz",compress = "gz")
ColonData$AllTypes <- ifelse(ColonData$RNA_snn_res.0.8 %in% 13,"MyeloidCells",ColonData$AllTypes)
ColonData$AllTypes <- ifelse(ColonData$RNA_snn_res.0.8 %in% c(19,24),"TCells",ColonData$AllTypes)
for(i in unique(ColonData$AllTypes)){
sub <- subset(ColonData,cells=rownames(ColonData@meta.data[ColonData@meta.data$AllTypes %in% i,]))
readr::write_rds(sub,path=paste("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/",i,"_harmony.rds.gz",sep=""),compress = "gz")
}
####Cluster distribution
EachSampleClusterDis <- lapply(unique(ColonDataS@meta.data$Pat_Tissues),function(Sample){
seur <- subset(ColonDataS,cells=rownames(ColonDataS@meta.data[ColonDataS@meta.data$Pat_Tissues %in% Sample,]))
seur@meta.data$AllTypes %>% table %>% data.frame %>% set_colnames(c("CellTypes","Number")) %>% dplyr::mutate(Per=100*Number/sum(Number))
})
names(EachSampleClusterDis) <- unique(ColonDataS@meta.data$Pat_Tissues)
EachSampleClusterDisA <- dplyr::bind_rows(EachSampleClusterDis) %>% dplyr::mutate(Sample=rep(names(EachSampleClusterDis),times=unlist(lapply(EachSampleClusterDis,nrow))))
pdf(paste("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/Cluster_Distribution.pdf",sep=""),width = 6,height = 4)
ggplot(data = EachSampleClusterDisA, aes(x = Sample, fill=CellTypes,y=Per)) + #log2(as.numeric(CopyNumber)))
geom_bar(stat = "identity",width = 0.6)+
scale_y_continuous(expand=c(0,0))+labs(y="percentage(%)")+
scale_x_discrete(limits=Samples)+
scale_fill_manual(values=.cluster_cols[c(9,1,12,3,8)],name="")+
theme(axis.text=element_text(color = "black"),
panel.background = element_blank(),panel.grid=element_blank(),
legend.title=element_blank(), axis.text.y=element_text(color = "black"),
axis.title.x = element_blank(),axis.line = element_line(color="black"),legend.text = element_text(size=10))
dev.off()
EachSampleClusterDisA[,c("Patients","Tissue")] <- data.frame(do.call(rbind,strsplit(as.character(EachSampleClusterDisA$Sample),"_")))
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/PercentageDefinedDistribution_Line.pdf",
width = 7,height = 3)
ggplot(data = EachSampleClusterDisA, aes(x = Tissue, color=Patients,y=Per,group=Patients)) + #log2(as.numeric(CopyNumber)))
geom_point(size=2)+
facet_wrap(~CellTypes,nrow=1,scale="free")+
scale_x_discrete(limit=c("N","T"),labels=c("N","T"))+
scale_y_continuous(expand=c(0.2,0))+
scale_color_manual(limit=c("p7","p9","p10","p11","p12"),values=RColorBrewer::brewer.pal(5,"Set2"),name="")+
geom_line(linetype = "dashed")+#stat_compare_means(label = "p.format")+
ylab("Percentage (%)")+
theme(axis.text=element_text(color = "black"),
panel.background = element_blank(),panel.grid=element_blank(),
legend.title=element_blank(), axis.text.y=element_text(color = "black"),
strip.background = element_blank(),strip.text = element_text(size=12),
axis.title.x = element_blank(),axis.line = element_line(color="black"),legend.text = element_text(size=10))+
ggpubr::stat_compare_means(comparisons = list(v1=c("N","T")),label="p.format",method = "t.test",paired =T )
dev.off()
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/Cluster_Distribution_Number.pdf",width = 8,height = 4)
ggplot(data = EachSampleClusterDisA, aes(x = Sample, y= CellTypes)) + #log2(as.numeric(CopyNumber)))
geom_tile(aes(fill=CellTypes),color="black")+
geom_text(aes(label=Number))+
coord_fixed(ratio = 0.8)+
scale_y_discrete(limits = c("MyeloidCells","TCells", "BCells","EpithelialCells", "StromalCells"))+labs(y="percentage(%)")+
scale_x_discrete(limits=names(EachSampleClusterDis))+
scale_fill_manual(values=.cluster_cols[c(9,1,12,3,8)],name="",guide=F)+
theme(axis.text=element_text(color = "black"),
panel.background = element_blank(),panel.grid=element_blank(),
axis.text.y=element_text(color = "black"),axis.ticks = element_blank(),
axis.title = element_blank(),legend.text = element_text(size=10))
dev.off()
###
####
#plasma cells CD38+; SDC1:CD138+
FeaturePlot(ColonData,reduction = "umap",features = c("CD19"),cols = c("gray","orange","red"))
2.majortypesfeaturedifference_all.r
library(magrittr)
library(ggplot2)
library(nlme)
library(Seurat)
options(stringsAsFactors = F)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/5.FunctionEnrichment/MarjorTypes/")
###
###
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineSubtype_AllFunctionTF.rds.gz")
###HallMark
Features <- c("RNA","Metab","BP","React", "Hallmark")
AllFeatureDiff <- lapply(Features,function(AssayName){
DefaultAssay(WholeTissueSC) <- AssayName
Idents(WholeTissueSC) <- "SubTypes"
if(AssayName=="RNA"){
DEfeatures <- FindAllMarkers(WholeTissueSC,min.pct=0.1,logfc.threshold = 0.25,pseudocount.use = 0.1,only.pos = T)
}else{
DEfeatures <- FindAllMarkers(WholeTissueSC,min.pct=0.1,logfc.threshold = 0.25,pseudocount.use = F,only.pos = T)
}
return(DEfeatures)
})
names(AllFeatureDiff ) <- Features
readr::write_rds(AllFeatureDiff,
path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/5.FunctionEnrichment/MarjorTypes/FeaturesDifferenceAmongSubtypes.rds.gz",compress = "gz")
SubTypesFeaturesDifference <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/5.FunctionEnrichment/MarjorTypes/SubTypesFeaturesDifference.rds.gz")
AllFeatureDiff <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/5.FunctionEnrichment/MarjorTypes/FeaturesDifferenceAmongSubtypes.rds.gz")
openxlsx::write.xlsx(AllFeatureDiff$RNA[AllFeatureDiff$RNA$avg_logFC > 1,],
file="DEgenesAmongSubtypes.xlsx")
2.subset_myeloid.r
library(Seurat)
library(ggplot2)
library(magrittr)
library(Matrix)
library(wordspace)
library(Hmisc)
library(rsvd)
options(stringsAsFactors = F)
set.seed(42)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/Myeloid/")
Samples <- c("p7_N","p9_N", "p10_N", "p11_N", "p12_N", "p7_T", "p9_T", "p10_T", "p11_T", "p12_T")
###Features
##1. mast cells, Expression profiling of constitutive mast cells reveals a unique identity within the immune system
##Mast: KIT (CD117) ,IL1RL1, Mast cells and basophils: trojan horses of conventional lin- stem/progenitor cell isolates
##Mast activated marker: ENPP3
.cluster_cols <- c(
"#DC050C", "#FB8072", "#1965B0", "#7BAFDE", "#882E72",
"#B17BA6", "#FF7F00", "#FDB462", "#E7298A", "#E78AC3",
"#33A02C", "#B2DF8A", "#55A1B1", "#8DD3C7", "#A6761D",
"#E6AB02", "#7570B3", "#BEAED4", "#666666", "#999999",
"#aa8282", "#d4b7b7", "#8600bf", "#ba5ce3", "#808000",
"#aeae5c", "#1e90ff", "#00bfff", "#56ff0d", "#ffff00")
TissueMyeloid <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/MyeloidCells_harmony.rds.gz")
#Idents(TissueMyeloid) <- TissueMyeloid$Pat_Tissues
TissueMyeloid <- NormalizeData(TissueMyeloid)
TissueMyeloid <- FindVariableFeatures(object =TissueMyeloid, dispersion.function = LogVMR)
TissueMyeloid <- ScaleData(object = TissueMyeloid)
TissueMyeloid <- RunPCA(object = TissueMyeloid, npcs = 30, verbose = FALSE)
TissueMyeloid <- TissueMyeloid %>%
RunHarmony("Pat_Tissues", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(TissueMyeloid, 'harmony')
TissueMyeloid <- TissueMyeloid %>%
RunUMAP(reduction = "harmony", dims = 1:20) %>%
FindNeighbors(reduction = "harmony", dims = 1:20)
TissueMyeloid <- TissueMyeloid %>%
RunTSNE(reduction = "harmony", dims = 1:20)
for(i in seq(0.1,0.8,by=0.1)){
TissueMyeloid <- FindClusters(TissueMyeloid, resolution =i, algorithm = 1)
}
TissueMyeloid@meta.data <- TissueMyeloid@meta.data[,c(1:6,grep("RNA_snn",colnames(TissueMyeloid@meta.data)))]
png("1.Cluster with Different resolution.png",width=3200,height=1400,res=200)
DiffResolution_plot <- lapply(seq(0.1,0.8,by=0.1),function(i){
p <- DimPlot( TissueMyeloid, reduction = "umap",label=T,group.by = paste("RNA_snn_res.",i,sep=""),cols = .cluster_cols)+
labs(title=paste("Resolution = ",i,sep=""))
return(p)
})
p_res <- cowplot::plot_grid(plotlist = DiffResolution_plot,ncol=4)
print(p_res)
dev.off()
DimPlot( TissueMyeloid, reduction = "umap",label=T,group.by = paste("RNA_snn_res.",0.8,sep=""),cols = .cluster_cols)
FeaturePlot(TissueMyeloid,features = c("ENPP3","IL1RL1","KIT"))
FeaturePlot(tt,features = c("CLEC4C","CLEC9A"))
tt <- readr::read_rds("~/WholeTissue/myeloid/Myeloid_tmp")
DimPlot(TissueMyeloid,group.by = "RNA_snn_res.0.8")
ttSingleR_Main <- SingleR::SingleR(GetAssayData(tt@assays$RNA, slot = "data"),
ref=ref.se,
method = "cluster",
clusters =tt@meta.data$RNA_snn_res.0.8,
labels = list(bp.se$label.main, hpca.se$label.main))
ttSingleR_Result <- data.frame(HPCA =ttSingleR_Main@listData$orig.results$HPCA$first.labels,
BP =ttSingleR_Main@listData$orig.results$BP$first.labels,
clusters=sort(unique(tt@meta.data$RNA_snn_res.0.8)))
ref.se <- readr::read_rds("~/Project/Collaboration/2020ShenLei/Matrix/ProcessData/SingeR_ref.se.rds.gz")
bp.se <- ref.se$BP
hpca.se <- ref.se$HPCA
TissueMyeloid_SingleR_Main <- SingleR::SingleR(GetAssayData(TissueMyeloid@assays$RNA, slot = "data"),
ref=ref.se,
method = "cluster",
clusters =TissueMyeloid@meta.data$RNA_snn_res.0.8,
labels = list(bp.se$label.main, hpca.se$label.main))
TissueMyeloidSingleR_Result <- data.frame(HPCA =TissueMyeloid_SingleR_Main@listData$orig.results$HPCA$first.labels,
BP =TissueMyeloid_SingleR_Main@listData$orig.results$BP$first.labels,
clusters=sort(unique(TissueMyeloid@meta.data$RNA_snn_res.0.8)))
TissueMyeloid@meta.data$HPCA <- SingleR_Result$HPCA[match(TissueMyeloid@meta.data$RNA_snn_res.0.8,SingleR_Result$clusters)]
TissueMyeloid@meta.data$BP <- SingleR_Result$BP[match(TissueMyeloid@meta.data$RNA_snn_res.0.8,SingleR_Result$clusters)]
####
MyeloidFake <- subset(TissueMyeloid,cells=rownames(TissueMyeloid@meta.data[TissueMyeloid@meta.data$RNA_snn_res.0.8 %in% 8,]))
readr::write_rds(MyeloidFake,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/MyeloidFake_Treal.rds.gz",compress = "gz")
MyeloidUpdate <- subset(TissueMyeloid,cells=rownames(TissueMyeloid@meta.data[!(TissueMyeloid@meta.data$RNA_snn_res.0.8 %in% 8),]))
readr::write_rds(MyeloidUpdate,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/MyeloidReal.rds.gz",compress = "gz")
###Macrophage
MacrophageMarker <- c("CD14", "FCGR3A", "FCGR1A", "CD68", "TFRC", "CCR5")
###Mast: KIT (CD117) ,IL1RL1, Mast cells and basophils: trojan horses of conventional lin- stem/progenitor cell isolates
###Mast activated marker: ENPP3
###Neutrophil:FCGR3B (CD16),CXCL8
##macrophages "C1Q" ##citation: Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma
# cDC1: c("IRF8","ID2","BATF3","NFIL3")
# cDC2: c("IRF4","KLF4","ZEB2")
##monocytes: "VCAN",
##DCs: "CLEC10A","CD1C"
##macrophages "C1Q" ##citation: Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma
###M2: "MRC1","CD163"
FeaturePlot(TissueMyeloid,features = c("CD14"))
2.subset_robustness.r
library(multiROC)
library(Seurat)
library(randomForest)
library(ggplot2)
library(caTools)
library(magrittr)
options(expressions = 5e5)
options(stringsAsFactors = F)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/Robustness")
Terms <- c("Endo","Epithelial","Myeloid","MSC","T","B")#
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz")
WholeTissueList <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineEachSubtypeUpdataB.rds.gz")
Clusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
Terms <- c("Endo","Epithelial","Myeloid","MSC","T","B")#
plotp <- list()
for(term in Terms){
Tg <- WholeTissueList[[term]]
TgCluster <- Clusters[Clusters$SubTypes %in% term,]
TgCluster$Cluster <- 1:nrow(TgCluster )
Tg <- ScaleData(Tg)
Tg_ScaleData <- Tg@assays$RNA@scale.data %>% t %>% data.frame
Tg_ScaleData$Cluster <- TgCluster$Cluster[match(Tg$DefineTypes,TgCluster$DefineTypes)]
#2 Get train and test data
set.seed(42)
split <- sample.split(Tg_ScaleData$Cluster, SplitRatio = 0.5)
train_df <- subset(Tg_ScaleData, split == TRUE)
test_df <- subset(Tg_ScaleData, split == FALSE)
#3 random frest
rf_res <- randomForest::randomForest(factor(Cluster) ~., data = train_df, ntree = 100)
rf_pred <- predict(rf_res, test_df, type = 'prob')
rf_pred <- data.frame(rf_pred)
colnames(rf_pred) <- paste("Cluster.",gsub("X","",colnames(rf_pred)), " _pred_RF",sep="")
#4 Merge true labels and predicted values
true_label <- dummies::dummy(test_df$Cluster, sep = ".")
true_label <- data.frame(true_label)
colnames(true_label) <- colnames(true_label)
colnames(true_label) <- paste(colnames(true_label), "_true")
final_df <- cbind(true_label, rf_pred)
#5 multiROC and multiPR
roc_res <- multi_roc(final_df, force_diag=T) #This function calculates the Specificity, Sensitivity and AUC of multi-class classifications.
pr_res <- multi_pr(final_df, force_diag=T) #his function calculates the Precision, Recall and AUC of multi-class classifications.
#6 Plot
plot_roc_df <- plot_roc_data(roc_res)
plot_pr_df <- plot_pr_data(pr_res)
plot_roc_df <- plot_roc_df[plot_roc_df$Group %in% paste("Cluster.",unique(test_df$Cluster)," ",sep=""), ]
plot_roc_df_AUC <- roc_res$AUC$RF %>% unlist %>% signif(digits = 3)
names(plot_roc_df_AUC) <- gsub(" ","",names(plot_roc_df_AUC))
plot_roc_df$Group <- gsub(" ","",plot_roc_df$Group)
plot_roc_df_AUC1 <- data.frame( plot_roc_df_AUC )
plot_roc_df_AUC1$Cluster <- rownames( plot_roc_df_AUC1)
plot_roc_df_AUC1 <- plot_roc_df_AUC1[!(plot_roc_df_AUC1$Cluster %in% c("macro","micro")),] %>% dplyr::mutate(x=1:nrow(.),y=1:nrow(.))
p <- ggplot(plot_roc_df_AUC1, aes(x = x, y=y,color=Cluster)) +
geom_point()+
labs(x="False Positive Rate",y="True Positive Rate")+
scale_color_manual(limit=paste("Cluster.",TgCluster$Cluster,sep=""),
values=as.character(TgCluster$Color),
label= paste(TgCluster$DefineTypes," (",plot_roc_df_AUC[match(paste("Cluster.",TgCluster$Cluster,sep=""),names(plot_roc_df_AUC))],")",sep=""),
name="AUC")+
theme_bw() + theme(plot.title = element_text(hjust = 0.5),
legend.justification=c(1, 0), legend.position=c(.95, .05),
legend.background = element_rect(fill=NULL, size=0.5,
linetype="solid", colour ="black"))
plotp[[term]] <- p
}
pdf("ClusterRobustness.pdf",width = 12,height = 8)
p <- cowplot::plot_grid(plotlist = plotp,ncol = 3)
print(p)
dev.off()
2.subset_stromal.r
library(Seurat)
library(ggplot2)
library(magrittr)
library(Matrix)
library(wordspace)
library(Hmisc)
library(rsvd)
options(stringsAsFactors = F)
set.seed(42)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/Stromal/")
Samples <- c("p7_N","p9_N", "p10_N", "p11_N", "p12_N", "p7_T", "p9_T", "p10_T", "p11_T", "p12_T")
###Features
##1. mast cells, Expression profiling of constitutive mast cells reveals a unique identity within the immune system
##Mast: KIT (CD117) ,IL1RL1, Mast cells and basophils: trojan horses of conventional lin- stem/progenitor cell isolates
##Mast activated marker: ENPP3
.cluster_cols <- c(
"#DC050C", "#FB8072", "#1965B0", "#7BAFDE", "#882E72",
"#B17BA6", "#FF7F00", "#FDB462", "#E7298A", "#E78AC3",
"#33A02C", "#B2DF8A", "#55A1B1", "#8DD3C7", "#A6761D",
"#E6AB02", "#7570B3", "#BEAED4", "#666666", "#999999",
"#aa8282", "#d4b7b7", "#8600bf", "#ba5ce3", "#808000",
"#aeae5c", "#1e90ff", "#00bfff", "#56ff0d", "#ffff00")
TissueStromal <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/StromalCells_harmony.rds.gz")
TissueStromal <- NormalizeData(TissueStromal)
TissueStromal <- FindVariableFeatures(object =TissueStromal, mean.function = ExpMean, dispersion.function = LogVMR)
TissueStromal <- ScaleData(object = TissueStromal)
TissueStromal <- RunPCA(object = TissueStromal, npcs = 30, verbose = FALSE)
TissueStromal <- TissueStromal %>%
RunHarmony("orig.ident", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(TissueStromal, 'harmony')
TissueStromal <- TissueStromal %>%
RunTSNE(reduction = "harmony", dims = 1:20)
TissueStromal <- TissueStromal %>%
RunUMAP(reduction = "harmony", dims = 1:30) %>%
FindNeighbors(reduction = "umap", dims = 1:2)
for(i in seq(0.05,0.3,by=0.05)){
TissueStromal <- FindClusters(TissueStromal, resolution =i, algorithm = 1)
}
TissueStromal@meta.data <- TissueStromal@meta.data[,c(1:6,grep("RNA_snn",colnames(TissueStromal@meta.data)))]
png("1.Cluster2 with Different resolution.png",width=2400,height=1400,res=200)
DiffResolution_plot <- lapply(seq(0.05,0.3,by=0.05),function(i){
p <- DimPlot( TissueStromal, reduction = "umap",label=T,group.by = paste("RNA_snn_res.",i,sep=""),cols = .cluster_cols)+
labs(title=paste("Resolution = ",i,sep=""))
return(p)
})
p_res <- cowplot::plot_grid(plotlist = DiffResolution_plot,ncol=3)
print(p_res)
dev.off()
ref.se <- readr::read_rds("~/Project/Collaboration/2020ShenLei/Matrix/ProcessData/SingeR_ref.se.rds.gz")
bp.se <- ref.se$BP
hpca.se <- ref.se$HPCA
TissueStromal_SingleR_Main <- SingleR::SingleR(GetAssayData(TissueStromal@assays$RNA, slot = "data"),
ref=ref.se,
method = "cluster",
clusters =TissueStromal@meta.data$RNA_snn_res.0.05,
labels = list(bp.se$label.main, hpca.se$label.main))
TissueStromalSingleR_Result <- data.frame(HPCA =TissueStromal_SingleR_Main@listData$orig.results$HPCA$first.labels,
BP =TissueStromal_SingleR_Main@listData$orig.results$BP$first.labels,
clusters=sort(unique(TissueStromal@meta.data$RNA_snn_res.0.05)))
TissueStromal@meta.data$HPCA <- TissueStromalSingleR_Result$HPCA[match(TissueStromal@meta.data$RNA_snn_res.0.05,TissueStromalSingleR_Result$clusters)]
TissueStromal@meta.data$BP <- TissueStromalSingleR_Result$BP[match(TissueStromal@meta.data$RNA_snn_res.0.05,TissueStromalSingleR_Result$clusters)]
DimPlot(TissueStromal,reduction = "umap",group.by="HPCA",cols = .cluster_cols)
DimPlot(TissueStromal,reduction = "umap",group.by="BP",cols = .cluster_cols)
pdf("Feature_Count_mito.pdf",width = 6,height = 9)
VlnPlot(object = TissueStromal, features = c("nCount_RNA","nFeature_RNA","percent.mito"),group.by="RNA_snn_res.0.05",pt.size = 0, cols= .cluster_cols,ncol = 1)+NoLegend()
dev.off()
FeaturePlot(TissueStromal,reduction = "umap",features = c("CD3E","CAFs",""),cols = c("gray","orange","red"))
###CAFs
CAFs <- c("FAP", "THY1", "DCN", "COL1A1", "COL1A2", "COL6A1", "COL6A2", "COL6A3")
TissueStromal@meta.data$CAFs <- apply(TissueStromal@assays$RNA@data[rownames(TissueStromal@assays$RNA@data) %in% CAFs,],2,mean)
pdf("CAFs_and markers.pdf",width = 9,height = 9)
p<-FeaturePlot(TissueStromal,reduction = "umap",features = c("CAFs",CAFs ),cols = c("gray","orange","red"),combine = F)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend()+theme(axis.title = element_blank(),axis.text = element_blank(),
axis.ticks = element_blank(),axis.line = element_line(color="gray",linetype = "dashed"))
}
cowplot::plot_grid(plotlist = p,ncol = 3)
dev.off()
panCAFs <- c("VIM", "COL1A1", "COL3A1")
TissueStromal@meta.data$panCAFs <- apply(TissueStromal@assays$RNA@data[rownames(TissueStromal@assays$RNA@data) %in% panCAFs,],2,mean)
FeaturePlot(TissueStromal,reduction = "umap",features = "panCAFs",cols = c("gray","orange","red"))
pdf("PanCAFs markers.pdf",width = 9,height = 3)
p<- FeaturePlot(TissueStromal,reduction = "umap",features = panCAFs,cols = c("gray","orange","red"),combine = F)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend()+theme(axis.title = element_blank(),axis.text = element_blank(),
axis.ticks = element_blank(),axis.line = element_line(color="gray",linetype = "dashed"))
}
cowplot::plot_grid(plotlist = p,ncol = 3)
dev.off()
Adipocytes <- c("ADH1B", "DLAT", "GPD1", "LBP", "PLIN1", "PPP1R1A","PPP2R1B", "PTGER3", "CILP", "ADIPOQ", "COL5A3", "PNPLA2")
TissueStromal@meta.data$Adipocytes <- apply(TissueStromal@assays$RNA@data[rownames(TissueStromal@assays$RNA@data) %in% Adipocytes,],2,mean)
FeaturePlot(TissueStromal,reduction = "umap",features = "Adipocytes",cols = c("gray","orange","red"))
FeaturePlot(TissueStromal,reduction = "umap",features = "PPP1R1A",cols = c("gray","orange","red"))
### myofibroblasts(MFs) citation: James Kinchen et al, Structural Remodeling of the Human Colonic Mesenchyme in Inflammatory Bowel Disease
MFs <- c("MYH11","ACTG2","MYOCD")
TissueStromal@meta.data$MFs <- apply(TissueStromal@assays$RNA@data[rownames(TissueStromal@assays$RNA@data) %in%MFs,],2,mean)
FeaturePlot(TissueStromal,reduction = "umap",features = "MFs",cols = c("gray","orange","red"))
FeaturePlot(TissueStromal,reduction = "umap",features = MFs,cols = c("gray","orange","red"))
pdf("MF markers.pdf",width = 8,height = 3)
p<- FeaturePlot(TissueStromal,reduction = "umap",features = MFs,cols = c("gray","orange","red"),combine = F)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend()+theme(axis.title = element_blank(),axis.text = element_blank(),
axis.ticks = element_blank(),axis.line = element_line(color="gray",linetype = "dashed"))
}
cowplot::plot_grid(plotlist = p,ncol = 3)
dev.off()
##Glial S100b, citation: James Kinchen et al, Structural Remodeling of the Human Colonic Mesenchyme in Inflammatory Bowel Disease
Glial <- c("S100B")
pdf("glial_marker_c8.pdf",width = 3,height = 3)
FeaturePlot(TissueStromal,reduction = "umap",features = "S100B",cols = c("gray","orange","red"))+NoLegend()
dev.off()
Pericytes <- c("RGS5")
pdf("Pericytes marker_c8.pdf",width = 3,height = 3)
FeaturePlot(TissueStromal,reduction = "umap",features = "RGS5",cols = c("gray","orange","red"))+NoLegend()
dev.off()
EndothelialCells <- c("PECAM1","VWF","CDH5")
TissueStromal@meta.data$EndothelialCells <- apply(TissueStromal@assays$RNA@data[rownames(TissueStromal@assays$RNA@data) %in% EndothelialCells,],2,mean)
FeaturePlot(TissueStromal,reduction = "umap",features = "EndothelialCells",cols = c("gray","orange","red"))
##mv endothelial cells: PLVAP
## endothelial WWRN1+
FeaturePlot(TissueStromal,reduction = "umap",features = "PLVAP",cols = c("gray","orange","red"))
FeaturePlot(TissueStromal,reduction = "umap",features = "ACKR1",cols = c("gray","orange","red"))
FeaturePlot(TissueStromal,reduction = "umap",features = "MMRN1",cols = c("gray","orange","red"))
pdf("Endothelial markers.pdf",width = 9,height =6)
p<- FeaturePlot(TissueStromal,reduction = "umap",features = c("PECAM1","VWF","CDH5","PLVAP","ACKR1","MMRN1"),cols = c("gray","orange","red"),combine = F)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend()+theme(axis.title = element_blank(),axis.text = element_blank(),
axis.ticks = element_blank(),axis.line = element_line(color="gray",linetype = "dashed"))
}
p <-FeaturePlot(ColonData,reduction = "umap",features = EndothelialCells,cols = c("gray","orange","orange","red","red"),combine = FALSE)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend() + NoAxes()
}
3.merge_clusterdefine.r
library(Seurat)
library(ggplot2)
library(magrittr)
library(harmony)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue")
SubTypesClusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/SubTypesColorCode.xlsx")
Clusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
TissuePlasmaB <- readr::read_rds("~/WholeTissue/BCell/PlasmaB/4thTissuePlasmacells.rds.gz")
TissuePlasmaB@meta.data <- TissuePlasmaB@meta.data[,c(1:6,ncol(TissuePlasmaB@meta.data))]
TissueB <- readr::read_rds("~/WholeTissue/BCell/BCellOnly/2ndTissueBcells.rds.gz")
TissueB@meta.data <- TissueB@meta.data[,c(1:6,ncol(TissueB@meta.data))]
ListB <- list(TissueB,TissuePlasmaB)
names(ListB) <- c("B","PlasmaB")
WholeTissueList <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineEachSubtypes.rds.gz")
WholeTissueList$Mast@meta.data <- WholeTissueList$Mast@meta.data[,1:7]
WholeTissueListNew <- c(WholeTissueList[c("Epithelial","MSC","Endo","Glia")],ListB ,WholeTissueList[c("T","ProliferatingT","ILCs","Myeloid","Mast")])
readr::write_rds(WholeTissueListNew,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineEachSubtypeUpdataB.rds.gz",compress = "gz")
openxlsx::write.xlsx(data.frame(DefineTypes=names(table(WholeTissueListNew$B$DefineTypes)),Colors=scales::hue_pal()(7)),file="tmp.xlsx")
ClustersColor <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
WholeTissueSC <- merge(x = WholeTissueListNew[[1]], y = WholeTissueListNew[names(WholeTissueListNew)[2:length(WholeTissueListNew)]],add.cell.ids=c(names(WholeTissueListNew)),merge.data=T,project = "SeuratProject")
ClusterNames <- lapply(WholeTissueListNew,function(x)names(table(x@meta.data$DefineTypes)))
MainTypes <- data.frame(Subtypes=names(WholeTissueListNew),MainTypes=c("Epithelial","Stroma","Stroma","Stroma","B","B","T","T","T","Myeloid","Myeloid"))
Clusters <- data.frame(DefineTypes=unlist(ClusterNames),SubTypes=rep(names(ClusterNames),times=unlist(lapply(ClusterNames,length))))
Clusters$MainTypes <- MainTypes$MainTypes[match(Clusters$SubTypes,MainTypes$Subtypes)]
WholeTissueSC@meta.data$SubTypes <- Clusters$SubTypes[match(WholeTissueSC@meta.data$DefineTypes,Clusters$DefineTypes)]
WholeTissueSC@meta.data$MainTypes <- MainTypes$MainTypes[match(WholeTissueSC@meta.data$SubTypes ,MainTypes$Subtypes)]
WholeTissueSC@meta.data$TopTypes <- ifelse(WholeTissueSC@meta.data$MainTypes %in% c("B","T","Myeloid"),"Immune",as.character(WholeTissueSC$MainTypes))
WholeTissueSC@meta.data$DefineTypes <- factor(WholeTissueSC@meta.data$DefineTypes,levels=Clusters$DefineTypes)
WholeTissueSC@meta.data$SubTypes <- ifelse(WholeTissueSC@meta.data$SubTypes %in% c("ILCs","ProliferatingT"),"T",as.character(WholeTissueSC@meta.data$SubTypes))
#WholeTissueSC@meta.data$SubTypes <- ifelse(WholeTissueSC@meta.data$DefineTypes %in% "Plasma","Plasma",as.character(WholeTissueSC@meta.data$SubTypes))
WholeTissueSC@meta.data$SubTypes <- factor(WholeTissueSC@meta.data$SubTypes,levels=c("Epithelial","MSC","Endo","Glia","B","PlasmaB","T","Myeloid","Mast"))
WholeTissueSC@meta.data$MainTypes <- factor(WholeTissueSC@meta.data$MainTypes,levels=c("Epithelial","Stroma","B","T","Myeloid"))
WholeTissueSC@meta.data$TopTypes <- factor(WholeTissueSC@meta.data$TopTypes,levels=c("Epithelial","Stroma","Immune"))
Patients <- data.frame(Pat_T = c("p1_N","p2_N", "p3_N", "p4_N", "p5_N", "p1_T", "p2_T", "p3_T", "p4_T", "p5_T"),
Pat_Tissues = c("p7_N","p9_N", "p10_N", "p11_N", "p12_N", "p7_T", "p9_T", "p10_T", "p11_T", "p12_T"))
WholeTissueSC <- NormalizeData(WholeTissueSC)
WholeTissueSC <- FindVariableFeatures(object =WholeTissueSC, mean.function = ExpMean, dispersion.function = LogVMR)
WholeTissueSC <- ScaleData(object = WholeTissueSC)
WholeTissueSC <- RunPCA(object = WholeTissueSC, npcs = 30, verbose = FALSE)
Idents(WholeTissueSC) <- WholeTissueSC$Pat_Tissues
WholeTissueSC <- WholeTissueSC %>%
RunHarmony("Pat_Tissues", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(WholeTissueSC, 'harmony')
WholeTissueSC <- WholeTissueSC %>%
RunUMAP(reduction = "harmony", dims = 1:20) %>%
FindNeighbors(reduction = "harmony", dims = 1:20)
WholeTissueSC <- WholeTissueSC %>%
RunTSNE(reduction = "harmony", dims = 1:10)
readr::write_rds(WholeTissueSC,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz",compress = "gz")
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz")
WholeTissueSC@meta.data$SubTypes <- factor(WholeTissueSC@meta.data$SubTypes,levels=SubTypesClusters$SubTypes)
tiff("WholeTissue_MainTypesDefineUpdateB.tiff",width = 1000,height = 1000,res=280)
DimPlot(WholeTissueSC,reduction = "umap",group.by = "SubTypes",cols = SubTypesClusters$Colors,pt.size = 0.1)+NoLegend()
dev.off()
pdf("ColorDefine.pdf",width = 4,height = 10)
ggplot(Clusters,aes(x=1,y=DefineTypes))+
geom_tile(aes(fill=DefineTypes))+
scale_y_discrete(limit=rev(Clusters$DefineTypes),labels=rev(Clusters$DefineTypes))+
scale_fill_manual(limit=Clusters$DefineTypes,values=Clusters$Colors,guide=F)+
coord_fixed()+theme(legend.background = element_blank(),panel.background = element_blank(),axis.text.x = element_blank(),
axis.title =element_blank(),axis.ticks = element_blank() )
dev.off()
pdf("ColorDefinePoint.pdf",width = 4,height = 10)
ggplot(Clusters,aes(x=1,y=DefineTypes))+
geom_point(aes(color=DefineTypes))+
scale_y_discrete(limit=rev(Clusters$DefineTypes),labels=rev(Clusters$DefineTypes))+
scale_color_manual(limit=Clusters$DefineTypes,values=Clusters$Colors,guide=F)+
coord_fixed()+theme(legend.background = element_blank(),panel.background = element_blank(),axis.text.x = element_blank(),
axis.title =element_blank(),axis.ticks = element_blank() )
dev.off()
pdf("SubTypesColorDefine1.pdf",width = 4,height = 4)
ggplot(SubTypesClusters,aes(x=1,y=SubTypes))+
geom_point(aes(color=SubTypes))+
scale_y_discrete(limit=rev(SubTypesClusters$SubTypes),labels=rev(SubTypesClusters$SubTypes))+
scale_color_manual(limit=SubTypesClusters$SubTypes,values=SubTypesClusters$Colors,guide=F)+
coord_fixed()+theme(legend.background = element_blank(),panel.background = element_blank(),axis.text.x = element_blank(),
axis.title =element_blank(),axis.ticks = element_blank() )
dev.off()
####
Idents(WholeTissueSC) <- WholeTissueSC$DefineTypes
WholeTissueSC_ClusterMarker <- FindAllMarkers(object = WholeTissueSC, only.pos = TRUE, min.pct = 0.2, logfc.threshold = 0.25)
openxlsx::write.xlsx(WholeTissueSC_ClusterMarker,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
WholeTissueSC_ClusterMarker <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
WholeTissueSC_ClusterMarkerF <- WholeTissueSC_ClusterMarker[WholeTissueSC_ClusterMarker$avg_logFC > 1,]
Idents(WholeTissueSC) <- WholeTissueSC$DefineTypes
WholeTissueSC_ClusterMarkerMast <- FindAllMarkers(object = WholeTissueSC, min.pct=0.1,logfc.threshold = 0.25,pseudocount.use = 0.1,only.pos = T,test.use = "MAST")
openxlsx::write.xlsx(WholeTissueSC_ClusterMarkerMast,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerTestByMAST.xlsx")
#####Plot Figures
WholeTissueListNew <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineEachSubtypeUpdataB.rds.gz")
readr::write_rds(WholeTissueListNew$Myeloid,path="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/CRC_Myeloid.rds.gz",compress = "gz")
p1 <- FeaturePlot(WholeTissueListA$MSC,feature="PLAU",split.by = "Tissues")
p2 <- FeaturePlot(WholeTissueSC,feature="ICAM1",split.by = "Tissues")
cowplot::plot_grid(plotlist = list(p1,p2))
3.merge_clusterdefine_plotfigure.r
library(Seurat)
library(ggplot2)
library(magrittr)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue")
PatientColor <- data.frame(Patients = c("p1","p2","p3","p4","p5"),PatColor=c("#D53E4F", "purple2", "#FDAE61", "#66C2A5", "#3288BD"))
##LRRC15
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz")
WholeTissueSC_ClusterMarker <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
WholeTissueSC_ClusterMarkerFC0.5 <- WholeTissueSC_ClusterMarker[WholeTissueSC_ClusterMarker$avg_logFC > 0.5,]
WholeTissueSC_ClusterMarkerFC0.5$DefineTypes <- WholeTissueSC_ClusterMarkerFC0.5$cluster
Clusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
SubClusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/SubTypesColorCode.xlsx")
SubClusters[,"SubTypes"] <- ifelse(SubClusters$SubTypes %in% "Plasma","PlasmaB",SubClusters$SubTypes)
DefineMarkers <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/AllCellTypeMarkers-0513.xlsx",sheet="DefineTypesMarkers")
PanMarkers <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/AllCellTypeMarkers-0513.xlsx",sheet="PanMarkers")
AddGenes <- setdiff(c(DefineMarkers$GeneSymbol,PanMarkers$GeneSymbol),WholeTissueSC_ClusterMarkerFC0.5$gene)
AddGenes <- DefineMarkers[DefineMarkers$GeneSymbol %in% AddGenes,] %>% dplyr::select(c(DefineTypes,GeneSymbol)) %>% set_colnames(c("DefineTypes","gene")) %>%
dplyr::mutate(avg_logFC=rep(1,times=nrow(.)))
DefineMarkers$DefineTypes <- factor(DefineMarkers$DefineTypes,levels=Clusters$DefineTypes)
PlotGenes <- unique(c(DefineMarkers$GeneSymbol,PanMarkers$GeneSymbol))
PlotGenesLeft <- unique(c(DefineMarkers[DefineMarkers$MainTypes %in% c("Epithelial","Stroma"),]$GeneSymbol,PanMarkers[PanMarkers$TopTypes %in% c("Epithelial","Stroma"),]$GeneSymbol))
PlotGenesRight <- unique(c(DefineMarkers[!(DefineMarkers$MainTypes %in% c("Epithelial","Stroma")),]$GeneSymbol,
PanMarkers[!(PanMarkers$TopTypes %in% c("Epithelial","Stroma")),]$GeneSymbol))
WholeTissueSC_ClusterMarkerFC0.5M <- WholeTissueSC_ClusterMarkerFC0.5 %>% dplyr::select(c(DefineTypes,gene,avg_logFC)) %>%
dplyr::bind_rows(AddGenes) %>% tidyr::spread(DefineTypes,avg_logFC,fill=0)
WholeTissueSC_ClusterMarkerFC0.5M <- WholeTissueSC_ClusterMarkerFC0.5M[,c("gene",Clusters$DefineTypes)]
#WholeTissueSC_ClusterMarkerFC0.5M <- WholeTissueSC_ClusterMarkerFC0.5M[!is.na(WholeTissueSC_ClusterMarkerFC0.5M$gene),Clusters$DefineTypes]
WholeTissueSC_ClusterMarkerFC0.5M <- WholeTissueSC_ClusterMarkerFC0.5M[do.call(order, c(WholeTissueSC_ClusterMarkerFC0.5M[,2:ncol(WholeTissueSC_ClusterMarkerFC0.5M)], list(decreasing=TRUE))),]
WholeTissue_ClusterID <- lapply(split(WholeTissueSC@meta.data,list(WholeTissueSC$DefineTypes)), function(x)rownames(x))
WholeTissueSCm <- WholeTissueSC@assays$RNA@data%>% as.matrix
colnames(WholeTissueSCm) <- gsub("\\.","-",colnames(WholeTissueSCm))
WholeTissueSCm <- WholeTissueSCm[apply(WholeTissueSCm,1,sum) > 0,]
#WholeTissueSCm <- WholeTissueSCm[rownames(WholetissueCom) %in% VariableFeatures(WholetissueHarmony),]
ClusterMean <- t(apply(WholeTissueSCm,1,function(x){
lapply(WholeTissue_ClusterID,function(ID){
mean(x[ID])
}) %>% unlist
}))
readr::write_rds(ClusterMean,path="EachClusterGeneMeanExpression.rds")
ClusterMean <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/EachClusterGeneMeanExpression.rds")
ClusterMean <- ClusterMean[apply(ClusterMean,1,sum)>0,]
ClusterMeanM <- ClusterMean[match(WholeTissueSC_ClusterMarkerFC0.5M$gene,rownames(ClusterMean)),]
ClusterMeanM[,1:ncol(ClusterMeanM)] <- t(apply(ClusterMeanM,1,scale))
Sample_order <- Clusters
ClusterMeanM[ClusterMeanM< -3] <- -3
ClusterMeanM[ClusterMeanM> 3] <- 3
Wholetissue_right_ha = ComplexHeatmap::rowAnnotation(foo = ComplexHeatmap::anno_mark(at = match(PlotGenes,WholeTissueSC_ClusterMarkerFC0.5M$gene), labels = PlotGenes ))
Wholetissue_RightPanel_ha = ComplexHeatmap::rowAnnotation(foo = ComplexHeatmap::anno_mark(at = match(PlotGenesRight,WholeTissueSC_ClusterMarkerFC0.5M$gene), labels = PlotGenesRight))
Wholetissue_LeftPanel_ha = ComplexHeatmap::rowAnnotation(foo = ComplexHeatmap::anno_mark(at = match(PlotGenesLeft,WholeTissueSC_ClusterMarkerFC0.5M$gene), labels = PlotGenesLeft))
Wholetissue_barcolor <- Clusters$Colors
names(Wholetissue_barcolor) <- Clusters$DefineTypes
Wholetissue_SampleBar <- Sample_order$DefineTypes
names(Wholetissue_SampleBar) <- Clusters$DefineTypes
Wholetissue_ha_bar <- ComplexHeatmap::HeatmapAnnotation(bar = Wholetissue_SampleBar,
col = list(bar = Wholetissue_barcolor))
pdf("Final_Wholetissue_HeatmapUpdateB.pdf",width = 8,height = 8)
ComplexHeatmap::Heatmap(ClusterMeanM,cluster_rows = FALSE,cluster_columns = FALSE,cluster_column_slices=F,use_raster = TRUE, name = "mat",
col=colorRampPalette(rev(RColorBrewer::brewer.pal(10,"RdYlBu"))[2:9],space="rgb")(100),#colorRampPalette(c("blue","white","red"),space="rgb")(100),#colorRampPalette(c("magenta", "black", "yellow"))(256),#
show_row_names = F,show_column_names = F,right_annotation = Wholetissue_right_ha,top_annotation = Wholetissue_ha_bar)
dev.off()
pdf("Final_Wholetissue_HeatmapUpdateB_Twoside.pdf",width = 8,height = 8)
ComplexHeatmap::Heatmap(ClusterMeanM,cluster_rows = FALSE,cluster_columns = FALSE,cluster_column_slices=F,column_order= names(Wholetissue_SampleBar),use_raster = TRUE, name = "mat",
col=colorRampPalette(rev(RColorBrewer::brewer.pal(10,"RdYlBu"))[2:9],space="rgb")(100),#colorRampPalette(c("blue","white","red"),space="rgb")(100),#colorRampPalette(c("magenta", "black", "yellow"))(256),#
show_row_names = F,show_column_names = F,
right_annotation = Wholetissue_RightPanel_ha,left_annotation =Wholetissue_LeftPanel_ha ,top_annotation = Wholetissue_ha_bar)
dev.off()
pdf("Final_Wholetissue_HeatmapMagentaYellow.pdf",width = 8,height = 8)
ComplexHeatmap::Heatmap(ClusterMeanM,cluster_rows = FALSE, cluster_columns = FALSE, cluster_column_slices=F,use_raster = TRUE, name = "mat",
col=colorRampPalette(c("magenta", "black", "yellow"))(256),#colorRampPalette(c("blue","white","red"),space="rgb")(100),#
show_row_names = F,show_column_names = F,right_annotation = Wholetissue_right_ha,top_annotation = Wholetissue_ha_bar)
dev.off()
ClustersColor <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
Tcolor <- ClustersColor[ClustersColor$SubTypes %in% "T",]
p1 <- DimPlot(WholeTissueList$T,reduction = "umap",group.by="DefineTypes",cols =Tcolor$Colors)+
coord_fixed()
p2 <- FeaturePlot(WholeTissueList$T,feature="GPR171",cols=c("gray93","blue"))+
coord_fixed()
p3 <- VlnPlot(WholeTissueList$T,feature="GPR171",group.by="DefineTypes",cols =Tcolor$Colors,pt.size = 0)+NoLegend()
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/GPR171_InCRC_Tcells.pdf",width = 6,height = 10)
cowplot::plot_grid(plotlist = list(p1,p2,p3),ncol=1,rel_heights = c(3.2,3.2,3.6))
dev.off()
SubTypesClusters <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/SubTypesColorCode.xlsx")
PP1 <- FeaturePlot(WholeTissueSC,feature="GPR171",cols=c("gray93","blue"))+
coord_fixed()
p1 <- DimPlot(WholeTissueSC,reduction = "umap",group.by="SubTypes",cols =SubTypesClusters$Colors)+
coord_fixed()
p3 <- VlnPlot(WholeTissueSC,feature="GPR171",group.by="SubTypes",cols =SubTypesClusters$Colors,pt.size = 0)+NoLegend()
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/GPR171_InCRC_Allcells.pdf",width = 5,height = 10)
cowplot::plot_grid(plotlist = list(p1,PP1,p3),ncol=1,rel_heights = c(3.2,3.2,3.6))
dev.off()
####
WholeTissueList <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineEachSubtypeUpdataB.rds.gz")
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz")
ClustersColor <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/ColorCode_Order.xlsx")
ClustersColor$MainTypes <- unique(WholeTissueSC@meta.data[,c("DefineTypes","MainTypes")])$MainTypes[match(ClustersColor$DefineTypes,unique(WholeTissueSC@meta.data[,c("DefineTypes","MainTypes")])$DefineTypes)]
#ClustersColor$DefineTypes <- factor(ClustersColor$DefineTypes,levels=names(table(WholeTissueSC$DefineTypes)))
#ClustersColor <- ClustersColor[order(ClustersColor$DefineTypes),]
DistributionPlot <- lapply(names(table(WholeTissueSC$MainTypes)),function(x){
SubSeurat <- subset(WholeTissueSC,cells=rownames(WholeTissueSC@meta.data[WholeTissueSC$MainTypes %in% x,]))
MSCcolor <- ClustersColor[ClustersColor$MainTypes %in% x,]
###All samples distribution.
SubSeurat@meta.data$DefineTypes <- factor(SubSeurat@meta.data$DefineTypes,levels=names(table(WholeTissueSC@meta.data[WholeTissueSC$MainTypes %in% x,]$DefineTypes)[table(WholeTissueSC@meta.data[WholeTissueSC$MainTypes %in% x,]$DefineTypes)>0]))
SubCellsDisA <- table(SubSeurat@meta.data[,c("Pat_Tissues","DefineTypes")]) %>%
data.frame %>% set_colnames(c("Sample","CellTypes","Number"))
SubCellsDisA <- lapply(split(SubCellsDisA,SubCellsDisA$Sample),function(X){
X%>% dplyr::mutate(Per=100*Number/sum(Number))
}) %>% dplyr::bind_rows(.)
SubCellsDisA[,c("PatientID","Tissues")] <- data.frame(do.call(rbind,strsplit(as.character(SubCellsDisA$Sample),"_")))
SubCellsDisA$Tissues <- factor(SubCellsDisA$Tissues,levels=c("N","T"),labels=c("Normal","Tumor"))
p <- ggplot(data = SubCellsDisA, aes(x = PatientID, fill=CellTypes,y=Per)) + #log2(as.numeric(CopyNumber)))
geom_bar(stat = "identity",width = 0.6)+
scale_y_continuous(expand=c(0,0))+labs(y="percentage(%)")+
facet_wrap(~Tissues)+
scale_x_discrete(limits=c(paste("p",1:5,sep="")))+
scale_fill_manual(limit=MSCcolor$DefineTypes,
values= MSCcolor$Colors,name="",guide=F)+
theme(axis.text.x=element_text(color = "black",hjust=0.5,size=9),
panel.background = element_blank(),panel.grid=element_blank(),strip.background = element_blank(),strip.text = element_text(size=12),
legend.title=element_blank(), axis.text.y=element_text(color = "black"),
axis.title.x = element_blank(),axis.line = element_line(color="black"),legend.text = element_text(size=8))
return(p)
})
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/DistributionPlot.pdf",width = 4,height = 20)
cowplot::plot_grid(plotlist = DistributionPlot,ncol=1)
dev.off()
Idents(WholeTissueSC) <- "SubTypes"
DistributionPlotData <- lapply(names(table(WholeTissueSC$SubTypes))[c(1:3,5,7:9)],function(x){
grep("Plasma",Cells(WholeTissueSC),value=T)
SubSeurat <- subset(WholeTissueSC,idents= x) #grep("Plasma",colnames(WholeTissueSC),value=T)
MSCcolor <- ClustersColor[ClustersColor$MainTypes %in% x,]
###All samples distribution.
SubSeurat@meta.data$DefineTypes <- factor(SubSeurat@meta.data$DefineTypes,
levels=names(table(WholeTissueSC@meta.data[WholeTissueSC$SubTypes %in% x,]$DefineTypes)[table(WholeTissueSC@meta.data[WholeTissueSC$SubTypes %in% x,]$DefineTypes)>0]))
SubCellsDisA <- table(SubSeurat@meta.data[,c("Pat_Tissues","DefineTypes")]) %>%
data.frame %>% set_colnames(c("Sample","CellTypes","Number"))
SubCellsDisA <- lapply(split(SubCellsDisA,SubCellsDisA$Sample),function(X){
X%>% dplyr::mutate(Per=100*Number/sum(Number))
}) %>% dplyr::bind_rows(.)
SubCellsDisA[,c("PatientID","Tissues")] <- data.frame(do.call(rbind,strsplit(as.character(SubCellsDisA$Sample),"_")))
SubCellsDisA$Tissues <- factor(SubCellsDisA$Tissues,levels=c("N","T"),labels=c("Normal","Tumor"))
return(SubCellsDisA)
})
names(DistributionPlotData) <- names(table(WholeTissueSC$SubTypes))[c(1:3,5,7:9)]
DiffCellTypes <- lapply(DistributionPlotData ,function(sub){
sub <- sub %>% dplyr::select(-c(Sample,Number) ) %>% tidyr::spread(Tissues,Per)
sapply(split(sub,sub$CellTypes),function(x){
aa <- data.frame(TumorMed=mean(x$Tumor),NormalMed=mean(x$Normal),Diff=mean(x$Tumor)-mean(x$Normal),
pval=t.test(x$Tumor,x$Normal,paired = T)$p.value)
return(aa)
}) %>% t %>% data.frame %>% dplyr::mutate(DefineTypes=rownames(.))
}) %>% dplyr::bind_rows()
DiffCellTypes$Diff <- as.numeric(DiffCellTypes$Diff)
DiffCellTypes <- merge(DiffCellTypes,unique(WholeTissueSC@meta.data[,c("SubTypes","DefineTypes")]),by="DefineTypes")
DiffCellTypes <- DiffCellTypes[order(DiffCellTypes$SubTypes,DiffCellTypes$Diff),]
openxlsx::write.xlsx(DiffCellTypes,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_InfiltrationDiff.xlsx")
DiffCellTypesSig <- DiffCellTypes[DiffCellTypes$Diff >0 & DiffCellTypes$pval < 0.05,]
WholeTissueSC_ClusterMarker <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
TumorUpMarkers <- WholeTissueSC_ClusterMarker[WholeTissueSC_ClusterMarker$avg_logFC > 0.75 & WholeTissueSC_ClusterMarker$cluster %in% DiffCellTypesSig$DefineTypes[grep("Prolifer",DiffCellTypesSig$DefineTypes,invert = T)],]
TumorUpMarkers$gene %>% unique
write.table(unique(TumorUpMarkers$gene),file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/TumorUpCellTypesMarkers.txt",quote = F,row.names = F,sep="\t")
DiffCellTypes$Class <- ifelse(DiffCellTypes$Diff >0 & DiffCellTypes$pval < 0.05,"Tumor",ifelse(DiffCellTypes$Diff < 0 & DiffCellTypes$pval < 0.05,"Normal","Both"))
DiffCellTypes$Diff <- signif(as.numeric(DiffCellTypes$Diff),digits = 3)
DiffCellTypes$Per <- -log2(as.numeric(DiffCellTypes$TumorMed)/as.numeric(DiffCellTypes$NormalMed))
DiffCellTypes$Class <- factor(DiffCellTypes$Class,levels=c("Normal","Both","Tumor"))
pdf("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/InfiltrationDifference.pdf",width = 4,height = 12)
ggplot(DiffCellTypes,aes(x=1,y=DefineTypes))+
geom_tile(aes(fill=Diff))+
scale_fill_gradientn(limit=c(-max(abs(DiffCellTypes$Diff)),max(abs(DiffCellTypes$Diff))),colours = colorRampPalette(colors=rev(RColorBrewer::brewer.pal(8,"RdYlBu")),space="rgb")(100))+
geom_text(aes(label=Diff))+
geom_tile(data=DiffCellTypes[DiffCellTypes$Class %in% "Tumor",],aes(x=1,y=DefineTypes),color="red",fill=NA)+
geom_tile(data=DiffCellTypes[DiffCellTypes$Class %in% "Normal",],aes(x=1,y=DefineTypes),color="blue",fill=NA)+
geom_tile(data=DiffCellTypes[DiffCellTypes$Class %in% "Both",],aes(x=1,y=DefineTypes),color="gray",fill=NA)+
scale_y_discrete(limit=DiffCellTypes[order(DiffCellTypes$Class,DiffCellTypes$Diff),]$DefineTypes)+
theme(axis.title = element_blank(),axis.ticks = element_blank(),panel.background = element_blank(),axis.text.x = element_blank())
dev.off()
openxlsx::write.xlsx(DiffCellTypes,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/DiffCellTypes_Normal_Tumor.xlxs")
###Patients colors
PatientUmapPlot <- lapply(c("Epithelial","MSC","Endo","T","Myeloid","Mast","B","PlasmaB"),function(x){
p <- DimPlot(WholeTissueList[[x]],reduction = "umap",group.by="orig.ident",cols = as.character(PatientColor$PatColor))+NoLegend()
return(p)
})
PatientUmapPlotN <- lapply(c("Epithelial","MSC","Endo","T","Myeloid","Mast","B","PlasmaB"),function(x){
sub <- WholeTissueList[[x]]
sub <- subset(sub,cells = rownames(sub@meta.data[sub@meta.data$Tissues %in% "N",]))
p <- DimPlot(sub,reduction = "umap",group.by="orig.ident",cols = as.character(PatientColor$PatColor))+NoLegend()
return(p)
})
PatientUmapPlotT <- lapply(c("Epithelial","MSC","Endo","T","Myeloid","Mast","B","PlasmaB"),function(x){
sub <- WholeTissueList[[x]]
sub <- subset(sub,cells = rownames(sub@meta.data[sub@meta.data$Tissues %in% "T",]))
p <- DimPlot(sub,reduction = "umap",group.by="orig.ident",cols = as.character(PatientColor$PatColor))+NoLegend()
return(p)
})
tiff("PatientUmapPlotSubTypes Plot.tiff",width = 2800,height = 1400,res=200)
cowplot::plot_grid(plotlist = PatientUmapPlot,ncol = 4)
dev.off()
tiff("PatientUmapPlotSubTypes PlotNT.tiff",width = 5600,height = 1400,res=200)
cowplot::plot_grid(plotlist = c(PatientUmapPlotT,PatientUmapPlotN),ncol = 8)
dev.off()
###Subtypes distribution
###MainTypes distribution
PatientUmapPlotNT <- lapply(c("MSC","Endo","T","Myeloid","Mast","B","PlasmaB","Epithelial"),function(x){
#sub <- subset(sub,cells = rownames(sub@meta.data[sub@meta.data$Tissues %in% "T",]))
p <- DimPlot(WholeTissueList[[x]],reduction = "umap",group.by="orig.ident",cols = as.character(PatientColor$PatColor),split.by = "Tissues")+NoLegend()
return(p)
})
names(PatientUmapPlotNT) <- c("MSC","Endo","T","Myeloid","Mast","B","PlasmaB","Epithelial")
tiff("PatientUmapPlotSubTypes PlotNT_.tiff",width = 1400,height = 4200,res=200)
cowplot::plot_grid(plotlist = PatientUmapPlotNT[c("Endo","T","Mast","B","PlasmaB","Epithelial")],ncol = 1)
dev.off()
tiff("PatientUmapPlotNT_.tiff",width = 800,height = 800,res=200)
DimPlot(WholeTissueSC,reduction = "umap",group.by="Tissues",cols = c("blue2","tomato"),pt.size = 0.1)+NoLegend()
dev.off()
WholeTissueList
WholeTissueSC_ClusterMarker <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
WholeTissueSC_ClusterMarkerF <- WholeTissueSC_ClusterMarker[WholeTissueSC_ClusterMarker$avg_logFC > 1,]
###Pan marker
PanMarkers <- c("EPCAM","COL1A1","COL3A1","CD19","MZB1","PTPRC","CD3E","CD14")
p <- FeaturePlot(WholeTissueSC,reduction = "umap",features = PanMarkers,cols = c("gray","orange","red"),pt.size = 0.5,combine = F)
for(i in 1:length(p)) {
p[[i]] <- p[[i]] + NoLegend()+theme(axis.title = element_blank(),axis.text = element_blank(),
axis.ticks = element_blank(),axis.line = element_line(color="gray",linetype = "dashed"))
}
tiff("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/PanMarker.tiff",width = 4200,height = 600,res=200)
cowplot::plot_grid(plotlist = p,ncol = 8)
dev.off()
4.0getmajortypesexpressionmatrix.r
library(corrplot)
library(survival)
library(ggpubr)
library(maxstat)
options(stringsAsFactors = F)
setwd("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/ImputeTCGA/CIBERSORTx/")
WholeTissueSC <- readr::read_rds("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/ProcessData/WholeTissue_DefineUpdateB.rds.gz")
WholeTissueSC_ClusterMarker <- openxlsx::read.xlsx("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/WholeTissue_SC_ClusterMarkerUpdateB.xlsx")
WholeTissueSC_ClusterMarkerFC0.5 <- WholeTissueSC_ClusterMarker[WholeTissueSC_ClusterMarker$avg_logFC > 0.5,]
NT_SCexpNorm <- WholeTissueSC@assays$RNA@data %>% as.matrix
T_SCexpNorm <- NT_SCexpNorm[,rownames(WholeTissueSC@meta.data[WholeTissueSC@meta.data$Tissues %in% "T" ,])]
T_SCexpNorm <- T_SCexpNorm[apply(T_SCexpNorm,1,sum) >1,]
T_SCexpNorm <- apply(T_SCexpNorm,2,function(x){signif(x,digits = 3)})
T_SCexpNorm <- T_SCexpNorm[rownames(T_SCexpNorm) %in% WholeTissueSC_ClusterMarker$gene,]
colnames(T_SCexpNorm) <- WholeTissueSC@meta.data$SubTypes[match(colnames(T_SCexpNorm),rownames(WholeTissueSC@meta.data))]
write.table(T_SCexpNorm,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/ImputeTCGA/Cibersortx_T_Major_Input.txt",quote = F,row.names = T,col.names = T,sep="\t")
#T_SCexpNorm <- read.delim("~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/ImputeTCGA/Cibersortx_Input_2020.txt")
NT_SCexpNorm <- WholeTissueSC@assays$RNA@data %>% as.matrix
N_SCexpNorm <- NT_SCexpNorm[,rownames(WholeTissueSC@meta.data[WholeTissueSC@meta.data$Tissues %in% "N" ,])]
N_SCexpNorm <- N_SCexpNorm[apply(N_SCexpNorm,1,sum) >1,]
N_SCexpNorm <- apply(N_SCexpNorm,2,function(x){signif(x,digits = 3)})
N_SCexpNorm <- N_SCexpNorm[rownames(N_SCexpNorm) %in% WholeTissueSC_ClusterMarker$gene,]
colnames(N_SCexpNorm) <- WholeTissueSC@meta.data$SubTypes[match(colnames(N_SCexpNorm),rownames(WholeTissueSC@meta.data))]
write.table(N_SCexpNorm,file="~/Project/Collaboration/2020ILC/WholeTissue/Harmony/Fig/3.WholeTissue/ImputeTCGA/Cibersortx_N_Major_Input.txt",quote = F,row.names = T,col.names = T,sep="\t")