# 导入Seurat包
library(Seurat)
library(ggplot2)
# 查看当前工作目录
getwd()

# 设置工作目录（将工作目录切换到指定路径下）
setwd("E:/单细胞示例")

# 读取10x数据，data.dir参数指定存放文件的路径
seurat_data <- Read10X(data.dir = "./GSE234527_RAW")

################################## 创建Seurat对象
{
  seurat_obj <- CreateSeuratObject(counts = seurat_data,
                                   project = "GSM7470392_352",
                                   min.features = 200,
                                   min.cells = 3)
  
  
  # 查看Seurat对象的基本信息
  seurat_obj
  
  ## Assays对象
  seurat_obj@assays
  seurat_obj@active.assay
  DefaultAssay(seurat_obj)
  ####查看meta.data信息
  meta <- seurat_obj@meta.data
  head(meta)
  
}
##################################数据质控
{
  ######添加线粒体信息
  seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")
  
  head(seurat_obj$percent.mt)
  ######添加核糖体基因信息
  seurat_obj[["percent.rb"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")
  
  head(seurat_obj$percent.rb)
  #####质控前可视化
  VlnPlot(seurat_obj, 
          fill.by = "feature", # "feature", "ident"
          features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"),
          ncol = 4, pt.size = 0.1) +
    theme(plot.title = element_text(size=10))
  
  #####相关性
  FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "percent.mt")
  FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
  FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "percent.rb")
  FeatureScatter(seurat_obj, feature1 = "percent.rb", feature2 = "percent.mt")
  
  
  #########过滤
  seurat_obj <- subset(seurat_obj, subset = 
                         percent.mt < 15 &  # 过滤线粒体基因占比>15%的细胞
                         nFeature_RNA > 500 &  # 过滤检测到基因数<500的细胞
                         nCount_RNA > 200)     # 过滤总UMI数<200的细胞
  VlnPlot(seurat_obj, 
          features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), 
          ncol = 4, pt.size = 0.1) +
    theme(plot.title = element_text(size=10))
  
}
