zl程序教程

您现在的位置是:首页 >  其它

当前栏目

失败经历的at2上皮细胞分群

失败 经历
2023-09-14 09:09:48 时间
library(Seurat)
library(ggplot2)
library(dplyr)
library(harmony)



load(file = "G:/silicosis/sicosis/yll/macrophage/AT2-cells/pure_as_inpaper/Krt+_Retnlg_ AT2 cells.rds")
#比例图
ggplot(subset_data@meta.data, aes(x=RNA_snn_res.0.06, fill=stim)) + geom_bar(position = "fill")
ggplot(subset_data@meta.data, aes(x=Idents(subset_data), fill=stim)) + geom_bar(position = "fill")
Krt_Retnlg_AT2_cells=subset_data


DotPlot(subset_data,features = c("Csf1r","Csf1","Csf2rb","Retnlg"))
FeaturePlot(subset_data,features = c("Csf1r","Csf1","Csf2rb"))
FeaturePlot(subset_data,features = c("Retnlg"))
DimPlot(subset_data,label = TRUE)


load("G:/silicosis/sicosis/yll/macrophage/no cluster2/0.3/pure_cluster3_in_allmerge-IM/silicosis_cluster_merge.rds")
table(All.merge$new.cluster.idents)
DimPlot(All.merge,label = T)
table(Idents(All.merge),All.merge$stim)


DimPlot(All.merge,label = TRUE)
DotPlot(All.merge,features = c("Igkv13-85","Ighv1-64","Hbb-bs","Jchain","Vim",
                               "Scgb3a1","Igha"))

#2 没有igha的所有
marcrophage_fibroblast=subset(All.merge,idents = c("IM","AM1",'AM2','AM3',
                                                   'Fibroblast'))

AT2_Macrophage_Fibroblast=merge(Krt_Retnlg_AT2_cells,marcrophage_fibroblast)
table(Idents(AT2_Macrophage_Fibroblast),AT2_Macrophage_Fibroblast$stim)


subset_data=AT2_Macrophage_Fibroblast
#重新分组 去除批次效应
if(1==1){#重新分组 去除批次效应
  subset_data[["percent.mt"]] <- PercentageFeatureSet(subset_data, pattern = "^mt-")
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data = subset_data %>%
    Seurat::NormalizeData(verbose = FALSE) %>%  
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
    ScaleData(verbose = FALSE) %>%
    RunPCA(npcs = 50, verbose = FALSE)
  ElbowPlot(subset_data, ndims = 50)
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
  harmony_embeddings <- Embeddings(subset_data, 'harmony') 
  #######################cluster
  dims = 1:30
  subset_data <- subset_data %>% 
    RunUMAP(reduction = "harmony", dims = dims) %>% 
    RunTSNE(reduction = "harmony", dims = dims) %>% 
    FindNeighbors(reduction = "harmony", dims = dims)
}
DimPlot(subset_data,label = TRUE)
getwd()
path="G:/silicosis/sicosis/yll/macrophage/AT2-cells/at2_macro_fibroblast_from_allmerge"
dir.create(path)
setwd(path)
#save(subset_data,file ="G:/silicosis/sicosis/yll/macrophage/AT2-cells/at2_macro_fibroblast_from_allmerge/at2_macro_fib_no_igha.rds" )
load("G:/silicosis/sicosis/yll/macrophage/AT2-cells/at2_macro_fibroblast_from_allmerge/at2_macro_fib_no_igha.rds")
at2_macro_fib_no_igha=subset_data

subset_data=at2_macro_fib_no_igha

#3
#所有的at2 和maco  fib
marcrophage_fibroblast_igha=subset(All.merge,idents = c("IM","AM1",'AM2','AM3',
                                                   'Fibroblast',
                                                   'Igha+ AT2  cell'))

all_at2_macro_fib=merge(marcrophage_fibroblast_igha,Krt_Retnlg_AT2_cells)
subset_data=all_at2_macro_fib
#重新分组 去除批次效应
if(1==1){#重新分组 去除批次效应
  subset_data[["percent.mt"]] <- PercentageFeatureSet(subset_data, pattern = "^mt-")
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data = subset_data %>%
    Seurat::NormalizeData(verbose = FALSE) %>%  
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
    ScaleData(verbose = FALSE) %>%
    RunPCA(npcs = 50, verbose = FALSE)
  ElbowPlot(subset_data, ndims = 50)
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
  harmony_embeddings <- Embeddings(subset_data, 'harmony') 
  #######################cluster
  dims = 1:30
  subset_data <- subset_data %>% 
    RunUMAP(reduction = "harmony", dims = dims) %>% 
    RunTSNE(reduction = "harmony", dims = dims) %>% 
    FindNeighbors(reduction = "harmony", dims = dims)
}

DimPlot(subset_data,label = TRUE)
subset_data$all_all_at2_macro_fib=Idents(subset_data)
getwd()
#save(subset_data,file = "G:/silicosis/sicosis/yll/macrophage/AT2-cells/at2_macro_fibroblast_from_allmerge/_all_at2_macro_fib.rds")
load("G:/silicosis/sicosis/yll/macrophage/AT2-cells/at2_macro_fibroblast_from_allmerge/_all_at2_macro_fib.rds")

FeaturePlot(subset_data,features = c("Igha",'Sftpc','Retnla','Retnlg',
                                     'Etv5', 'Lamp3', 'Slc34a2',
 'Bex2', 'Cd36', 'Chi3l1', 'Chsy1', 'Cxcl15', 'Dlk1', 'Egfl6', 'Fabp12', 'Fabp5',
 'Fasn', 'Glrx', 'Hc', 'Il33',  'Lcn2', 'Lgi3', 'Lyz1', 'Lyz2', 'Mid1ip1', 
 'Mlc1', 'Rab27a', 'Retnla', 'S100g', 'Scd1', 'Sftpa1', 'Sftpb', 'Sftpc', 'Soat1'))

getwd()
pdf("at.pdf",width = 9,height = 19)
print(p)
dev.off()

DotPlot(subset_data,features = c("Igha",'Sftpc','Retnla','Retnlg',
                                     'Etv5', 'Lamp3', 'Slc34a2',
                                     'Bex2', 'Cd36', 'Chi3l1', 'Chsy1', 'Cxcl15', 'Dlk1', 'Egfl6', 'Fabp12', 'Fabp5',
                                     'Fasn', 'Glrx', 'Hc', 'Il33',  'Lcn2', 'Lgi3', 'Lyz1', 'Lyz2', 'Mid1ip1', 
                                     'Mlc1', 'Rab27a', 'S100g', 'Scd1', 'Sftpa1', 'Sftpb', 'Soat1'))

DotPlot(subset_data,features = c("Krt8",'Krt18',"Retnlg","Retnla","Sftpc",'Lamp3'))
FeaturePlot(subset_data,features = c("Krt8",'Krt18',"Retnlg","Retnla",
                                     "Sftpc","Sftpa1",'Sftpb'))
table(Idents(subset_data),subset_data$stim)
DimPlot(subset_data,label = TRUE)
Seurat::DimPlot(subset_data,split.by = 'stim',label = TRUE)

subset_data=FindClusters(subset_data)
DimPlot(subset_data,label = TRUE)





#4 全部上皮
load("G:/silicosis/sicosis/yll/macrophage/no cluster2/0.3/pure_cluster3_in_allmerge-IM/silicosis_cluster_merge.rds")
table(All.merge$new.cluster.idents)
DimPlot(All.merge,label = TRUE)

endo_at2=subset(All.merge,idents = c('Endothelial cell-2','Endothelial cell-1',
                                                   'AT2 cell-1','AT2 cell-2','Igha+ AT2  cell'))


subset_data=endo_at2
#重新分组 去除批次效应
if(1==1){#重新分组 去除批次效应
  subset_data[["percent.mt"]] <- PercentageFeatureSet(subset_data, pattern = "^mt-")
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data = subset_data %>%
    Seurat::NormalizeData(verbose = FALSE) %>%  
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
    ScaleData(verbose = FALSE) %>%
    RunPCA(npcs = 50, verbose = FALSE)
  ElbowPlot(subset_data, ndims = 50)
  VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
  harmony_embeddings <- Embeddings(subset_data, 'harmony') 
  #######################cluster
  dims = 1:10
  subset_data <- subset_data %>% 
    RunUMAP(reduction = "harmony", dims = dims) %>% 
    RunTSNE(reduction = "harmony", dims = dims) %>% 
    FindNeighbors(reduction = "harmony", dims = dims)
}

DimPlot(subset_data,label = TRUE)
subset_data$endo_At2=Idents(subset_data)
getwd()

subset_data=FindClusters(subset_data)
DimPlot(subset_data,label = TRUE)


subset_data=subset(subset_data,idents = seq(1,17,1))
DimPlot(subset_data,label = TRUE)
at2markers=FindAllMarkers(subset_data,only.pos = TRUE,logfc.threshold = 0.3)