• cross-tissue 成纤维细胞比例.r


    
    getwd()
    setwd("G:\\lung fibrosis")
    
    ##########mouse stable state--------------
    mouse_ssfibro=readRDS("./Mouse_SS_Fibro.RDS")
    library(Seurat)
    head(mouse_ssfibro@meta.data)
    DimPlot(mouse_ssfibro,label = TRUE,raster=FALSE)
    table(mouse_ssfibro$Tissue)
    dim(mouse_ssfibro)
    
    Idents(mouse_ssfibro)=mouse_ssfibro$Tissue
    
    lung_ss=subset(mouse_ssfibro,idents = "Lung")
    DimPlot(lung_ss,label = TRUE)
    
    str(lung_ss)
    
    subset_data=lung_ss
    subset_data[["percent.mt"]] <- PercentageFeatureSet(subset_data, pattern = "^mt-")
    grep(pattern = "mt",x = rownames(subset_data),value = TRUE)
    
    VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
    
    library(dplyr)
    subset_data = subset_data %>%
      Seurat::NormalizeData(verbose = FALSE) %>%  
      FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
      ScaleData(verbose = FALSE) %>%
      RunPCA(npcs = 50, verbose = FALSE)
    
    
    dims = 1:30
    subset_data <- subset_data %>% 
      RunUMAP(reduction = "pca", dims = dims) %>% 
      RunTSNE(reduction = "pca", dims = dims) %>% 
      FindNeighbors(reduction = "pca", dims = dims)
    
    
    #subset_data@meta.data$stim <- c(rep("Exp", length(grep("1$",colnames(subset_data)))),rep("Con", length(grep("2$",colnames(subset_data)))))
    table(subset_data)
    DimPlot(subset_data,label = TRUE,label.size = 7)+
      ggtitle("stable_state_mouse_lung")
    
    
    DotPlot(subset_data,features = c("Inmt","Cd34","Sepp1","Ly6a",
                                     "Gpx3","Selenop"))
    library(stringr)
    library(ggplot2)
    DotPlot(subset_data,features = str_to_title(c("GPX1","GPX3",
                                            "Dio1","Gpx4","Rps14","Selh","Selk","Selm","Selo","Sels","Selt",
                                            "Sep15","Sepp1","Sepw1","Sepx1","Txnrd1",
                                            "Txnrd2","Txnrd3")))+
      RotatedAxis()+ggtitle("stable_state_mouse_lung")
    ##selected stable lung--------
    dim(subset_data)
    
    
    
    ################mouse ps fibro----------
    mouse_ps_fibro=readRDS("./Mouse_PS_Fibro.RDS")
    head(mouse_ps_fibro@meta.data)
    table(mouse_ps_fibro$Tissue)
    Idents(mouse_ps_fibro)=mouse_ps_fibro$Tissue
    
    lung_fibro=subset(mouse_ps_fibro,idents = c('Lung Bleo_Fibrosis D14',
                                            "Lung Fibrosis D21",  #"Lung Nin_Fibrosis D11",
                                            "Lung Bleo_Fibrosis D11"))
    DimPlot(lung_fibro,label = TRUE)
    
    mouse_selected_lung_ps=lung_fibro
    {
      mouse_selected_lung_ps[["percent.mt"]] <- PercentageFeatureSet(mouse_selected_lung_ps, pattern = "^mt-")
      VlnPlot(mouse_selected_lung_ps, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
      mouse_selected_lung_ps = mouse_selected_lung_ps %>%
        Seurat::NormalizeData(verbose = FALSE) %>%  
        FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
        ScaleData(verbose = FALSE) %>%
        RunPCA(npcs = 50, verbose = FALSE)
      table(Idents(mouse_selected_lung_ps))
      mouse_selected_lung_ps$stim=Idents(mouse_selected_lung_ps)
      #mouse_selected_lung_ps@meta.data$stim <- c(rep("Exp", length(grep("1$",colnames(mouse_selected_lung_ps)))),rep("Con", length(grep("2$",colnames(mouse_selected_lung_ps)))))
      table(mouse_selected_lung_ps$stim)
      
      ##########################run harmony
      #BiocManager::install('harmony')
      library('harmony')
      
      mouse_selected_lung_ps <- mouse_selected_lung_ps %>% RunHarmony("stim", plot_convergence = TRUE)
      harmony_embeddings <- Embeddings(mouse_selected_lung_ps, 'harmony') 
      #######################cluster
      dims = 1:30
      mouse_selected_lung_ps <- mouse_selected_lung_ps %>% 
        RunUMAP(reduction = "harmony", dims = dims) %>% 
        RunTSNE(reduction = "harmony", dims = dims) %>% 
        FindNeighbors(reduction = "harmony", dims = dims)%>% 
        FindClusters()
      
      
      DimPlot(mouse_selected_lung_ps,label = TRUE,label.size = 7)+
        ggtitle("mouse_selected_lung_ps")
      
    }
    
    DotPlot(mouse_selected_lung_ps,features = str_to_title(c("GPX1","GPX3",
                                                  "Dio1","Gpx4","Rps14","Selh","Selk","Selm","Selo","Sels","Selt",
                                                  "Sep15","Sepp1","Sepw1","Sepx1","Txnrd1",
                                                  "Txnrd2","Txnrd3")))+
      RotatedAxis()+ggtitle("mouse_selected_lung_ps")
    
    
    DotPlot(mouse_selected_lung_ps,features = c("Inmt","Cd34","Sepp1","Ly6a",
                                     "Gpx3","Selenop","Selp"))
    
    dim(mouse_selected_lung_ps)
    #######merge stable and perturbed- stable and fibrosi-------
    subset_data$stim="stable"
    
    All.merge=merge(subset_data,mouse_selected_lung_ps)
    
    head(All.merge@meta.data)
    table(All.merge$stim)
    table(is.na(All.merge$stim))
    All.merge$stim=case_when(is.na(All.merge$stim)~"stable",
                             .default = All.merge$stim)
    table(All.merge$stim)
    
    
    {
      
      All.merge[["percent.mt"]] <- PercentageFeatureSet(All.merge, pattern = "^mt-")
      VlnPlot(All.merge, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
      All.merge = All.merge %>%
        Seurat::NormalizeData(verbose = FALSE) %>%  
        FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
        ScaleData(verbose = FALSE) %>%
        RunPCA(npcs = 50, verbose = FALSE)
      
      #All.merge@meta.data$stim <- c(rep("Exp", length(grep("1$",colnames(All.merge)))),rep("Con", length(grep("2$",colnames(All.merge)))))
      table(All.merge$stim)
      
      library('harmony')
      
      All.merge <- All.merge %>% RunHarmony("stim", plot_convergence = TRUE)
      harmony_embeddings <- Embeddings(All.merge, 'harmony') 
      #######################cluster
      dims = 1:30
      All.merge <- All.merge %>% 
        RunUMAP(reduction = "harmony", dims = dims) %>% 
        RunTSNE(reduction = "harmony", dims = dims) %>% 
        FindNeighbors(reduction = "harmony", dims = dims)
      
      All.merge=FindClusters(All.merge)
      DimPlot(All.merge,label = TRUE,label.size = 7)
    }
    
    All.merge=All.merg
    library(stringr)
    DotPlot(All.merge,features = str_to_title(c("GPX1","GPX3",
                                                             "Dio1","Gpx4","Rps14","Selh","Selk","Selm","Selo","Sels","Selt",
                                                             "Sep15","Sepp1","Sepw1","Sepx1","Txnrd1",
                                                             "Txnrd2","Txnrd3")))+
      RotatedAxis()+ggtitle("allmerge_stable and fibrosis")
    
    
    DotPlot(All.merge,features = c("Inmt","Cd34","Sepp1","Ly6a",
                                                "Gpx3","Selenop","Selp"))
    
    DimPlot(All.merge,label = TRUE,
            split.by = "stim")
    table(All.merge$stim)
    All.merge$group =ifelse(All.merge$stim=="stable","Stable","Fibrosis")
    
    FeaturePlot(All.merge,features = c("Inmt","Cd34","Sepp1","Ly6a",
                                   "Gpx3","Selenop"),split.by = "group")
    
    
    grep(pattern = "sele|sep|sel|inmt|selen",rownames(All.merge),ignore.case = TRUE,
         value = TRUE)
    
    All.merg=RenameIdents(All.merge,"2"="universal fib",
                         "5"="universal fib",
                         "13"="universal fib",
                         "0"="specialized fib",
                         "3"="specialized fib")
    ggplot(All.merg@meta.data, 
           aes(x=Idents(All.merg),fill=Idents(All.merg))) + geom_bar(position = "fill")  
     #+scale_fill_manual(values = c("royalblue", "red"))#scale_fill_gradientn(colors = colorRampPalette(c("royalblue", "firebrick3"))(12))#colorRampPalette(c("royalblue","firebrick3"))(2)
    
    ggplot(All.merg@meta.data, 
           aes(x=Idents(All.merg), fill=group)) + geom_bar(position = "fill")+
      RotatedAxis()
    
      table(All.merg$group,Idents(All.merg))  #%>%summarise()
    table(All.merg$group)
    
    getwd()
    save(All.merg,file = "All.merg.rds")
    
    
    dir.create("G:/lung fibrosis/silicosis_cross_tissue_fibroblasts")
    setwd("G:/lung fibrosis/silicosis_cross_tissue_fibroblasts")
    load("G:/lung fibrosis/All.merg.rds")
    library(ggplot2)
    library(Seurat)
    head(All.merg@meta.data)
    table(All.merg$stim)
    table(All.merg$group)
    Idents(All.merg)=All.merg$group
    
    degs_cross_tissue_mice_fibrosis_vs_stable=FindMarkers(All.merg,
                                                          ident.1 = "Fibrosis",
                                                            ident.2 = "Stable")
    degs_cross_tissue_mice_fibrosis_vs_stable$gene=rownames(degs_cross_tissue_mice_fibrosis_vs_stable)
    
    DT::datatable(degs_cross_tissue_mice_fibrosis_vs_stable)
    openxlsx::write.xlsx(degs_cross_tissue_mice_fibrosis_vs_stable,
                         file ="degs_cross_tissue_mice_fibrosis_vs_stable.xlsx" )
    
    
    
    
    
    
    
    • 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
    • 151
    • 152
    • 153
    • 154
    • 155
    • 156
    • 157
    • 158
    • 159
    • 160
    • 161
    • 162
    • 163
    • 164
    • 165
    • 166
    • 167
    • 168
    • 169
    • 170
    • 171
    • 172
    • 173
    • 174
    • 175
    • 176
    • 177
    • 178
    • 179
    • 180
    • 181
    • 182
    • 183
    • 184
    • 185
    • 186
    • 187
    • 188
    • 189
    • 190
    • 191
    • 192
    • 193
    • 194
    • 195
    • 196
    • 197
    • 198
    • 199
    • 200
    • 201
    • 202
    • 203
    • 204
    • 205
    • 206
    • 207
    • 208
    • 209
    • 210
    • 211
    • 212
    • 213
    • 214
    • 215
    • 216
    • 217
    • 218
    • 219
    • 220
    • 221
    • 222
    • 223
    • 224
    • 225
  • 相关阅读:
    LeetCode136只出现一次的数字
    阿里巴巴最新的 SpringCloud 面试题是如何让我“颜面扫地”的?
    boost之跨平台 错误处理
    反弹shell原理与实现
    分布式调度 Elastic-job
    电脑不小心格式化了,怎么恢复?
    第一章《初学者问题大集合》第8节:Java程序初步释疑
    Spring web MVC(入门)
    从执行class文件开始认识JVM
    3月14日,每日信息差
  • 原文地址:https://blog.csdn.net/qq_52813185/article/details/132893178