基于单细胞和整体RNA测序的氧化应激反应生物标志物的卵巢癌研究（文章复现）3.单细胞分析

[ixxmu]

https://mp.weixin.qq.com/s/oi4M6UetEjepu83ux8jR9A

[ixxmu]

基于单细胞和整体RNA测序的氧化应激反应生物标志物的卵巢癌研究（文章复现）3.单细胞分析 by R 学习践行者

单细胞RNA测序数据集GSE146026的表达矩阵和元数据信息文件从（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE146026）数据库中获得。下载好数据放置在工作目录上。

单细胞测序数据分析流程：

      使用“Seurat” R 包对单细胞测序数据进行了分析。分析主要包括以下步骤：构建对象、数据标准化、数据降维聚类和识别标记基因。使用CreateSeuratObject构建Seurat对象。将最小细胞数设为3，最小特征数设为200，并首次对基因和细胞进行过滤。以下阈值用于再次过滤细胞：特征数为（200，6000），计数数为（300，40000），过表达double_cutoff值为0.15，线粒体基因比例< 30%，红细胞读取比例< 20%。选择前2000个高变基因作为主成分分析（PCA）的输入。通过ElbowPlot函数确定前15个主成分（PCs）为后续分析的重要主成分。

       由于单细胞测序数据来自多个数据集，因此使用“harmony” R 包进行批次校正，以消除干扰下游分析的批次效应。调试后并参考原始贡献中的聚类结果，使用FindClusters算法以分辨率0.5识别肿瘤细胞子集。

      统一流形近似与投影（UMAP）是一种基于黎曼几何生成数据理论框架的降维流行学习技术。UMAP在处理高通量和高维数据（如单细胞测序）时保留了更多的全局结构，性能优越，扩展性更好。使用UMAP对PCA获得的数据进一步降维。使用DimPlot函数将细胞类型划分为低维空间并进行可视化。使用FeaturePlot、DotPlot、DoHeatmap和VlnPlot函数可视化肿瘤细胞的表达、分布景观和特征基因。使用DotPlot函数可视化每个细胞亚群中与特定氧化应激反应相关的标记基因。选择|avg_log2FC| > 3的簇标记基因和氧化应激反应因子的交集基因中的顶级基因进行映射显示。

1.创建Seurat对象

library(tidyverse)
a0 = data.table::fread("GSE146026_Izar_HGSOC_ascites_10x_log.tsv.gz",
                      data.table = F) #读取数据

a = column_to_rownames(a0,var = "Cell_ID") #列变成行名
a[1:10,1:2]

##                                    10x_1                        10x_2
## 10x_barcode 10x_3288_t1_AAACATACCTTCCG-1 10x_3288_t1_AAACATACTCCTAT-1
## patient                                5                            5
## time                                   1                            1
## sample_ID                         3288.1                       3288.1
## clst                                   1                            1
## TSNE_x                      4.533631e+01                 3.507609e+01
## TSNE_y                      4.693348e+01                -2.010105e+01
## AL627309.1                      0.000000                     0.000000
## LINC00115                       0.000000                     0.000000
## SAMD11                         0.0000000                    0.0000000

前7行是基本信息，不是表达数据，所以需要整理一下。

exp0=a[-(1:7),]
exp0 = apply(exp0,2,as.numeric)
rownames(exp0) = rownames(a)[-(1:7)] 
exp0[1:4,1:4]

##            10x_1 10x_2 10x_3    10x_4
## AL627309.1     0     0     0 0.000000
## LINC00115      0     0     0 0.000000
## SAMD11         0     0     0 0.000000
## NOC2L          0     0     0 3.598553

#获取元数据的信息，比如细胞类型或样本来源
meta = as.data.frame(t(a[1:7,]))

library(Seurat)
sce.all <- CreateSeuratObject(
  counts = exp0,  #创建Seurat对象
  project = "ROS", 
  meta.data = meta,
  min.cells = 3,  #至少有3个细胞表达某个基因，才能保留这个基因
  min.features = 200) #每个细胞至少表达200个基因，才能保留这个细胞

#从Seurat对象中提取RNA计数矩阵
exp = sce.all[["RNA"]]@layers$counts
dim(exp)

## [1] 11548  9609

exp[1:4,1:4]

## 4 x 4 sparse Matrix of class "dgCMatrix"
##                    
## [1,] . . . .       
## [2,] . . . .       
## [3,] . . . .       
## [4,] . . . 3.598553

[email protected][1,] #元数据的信息

##       orig.ident nCount_RNA nFeature_RNA                  10x_barcode patient
## 10x_1        10x   4677.995          839 10x_3288_t1_AAACATACCTTCCG-1       5
##       time sample_ID clst       TSNE_x       TSNE_y
## 10x_1    1    3288.1    1 4.533631e+01 4.693348e+01

#将患者信息赋给orig.ident
sce.all$orig.ident = sce.all$patient
table(sce.all$orig.ident) #查看每个患者的细胞数量分布

## 
##    1    2    3    4    5    6 
##  354  573  821  671 5998 1192

2.质控

sce.all[["percent.mt"]] <- PercentageFeatureSet(sce.all, 
                                                pattern = "^MT-") #计算线粒体基因的百分比
sce.all[["percent.rp"]] <- PercentageFeatureSet(sce.all, 
                                                pattern = "^RP[SL]") #计算核糖体蛋白基因的百分比
sce.all[["percent.hb"]] <- PercentageFeatureSet(sce.all, 
                                                pattern = "^HB[^(P)]") #计算血红蛋白基因的百分比

head([email protected], 1) #元信息表会添加计算出来的百分数

##       orig.ident nCount_RNA nFeature_RNA                  10x_barcode patient
## 10x_1          5   4677.995          839 10x_3288_t1_AAACATACCTTCCG-1       5
##       time sample_ID clst       TSNE_x       TSNE_y percent.mt percent.rp
## 10x_1    1    3288.1    1 4.533631e+01 4.693348e+01   1.820066   7.208477
##       percent.hb
## 10x_1          0

使用VlnPlot函数绘制小提琴图，以可视化每个细胞样本的基因特征（nFeature_RNA）、RNA计数（nCount_RNA）、线粒体基因百分比（percent.mt）、核糖体蛋白基因百分比（percent.rp）和血红蛋白基因百分比（percent.hb）。

VlnPlot(sce.all, 
        features = c("nFeature_RNA",
                     "nCount_RNA", 
                     "percent.mt",
                     "percent.rp",
                     "percent.hb"), 
        ncol = 3, #3列
        pt.size = 0,  #数据点大小
        group.by = "orig.ident") #按orig.ident分组

#对单细胞RNA测序数据进行过滤，以保留符合特定条件的细胞
sce.all = subset(sce.all,
                percent.mt < 5 &
                nFeature_RNA < 6200 &
                nCount_RNA < 18000 &
                percent.rp <15 &
                percent.hb <1)

1）筛选线粒体基因百分比（percent.mt < 5）：

只保留线粒体基因表达百分比低于5%的细胞。高线粒体基因百分比通常表示细胞质量较低，可能是因为细胞正在经历应激或死亡。

2）筛选基因特征数（nFeature_RNA < 6200）：

 只保留基因特征数少于6200的细胞。基因特征数过高可能表示双重细胞（doublets）或其他技术噪音。

3）筛选RNA计数（nCount_RNA < 18000）：

 只保留RNA计数少于18000的细胞。RNA计数过高可能也表示双重细胞或技术噪音。

4）筛选核糖体蛋白基因百分比（percent.rp < 15）：

 只保留核糖体蛋白基因表达百分比低于15%的细胞。核糖体蛋白基因的高表达可能会影响下游分析的结果。

5）筛选血红蛋白基因百分比（percent.hb < 1）：

 只保留血红蛋白基因表达百分比低于1%的细胞。血红蛋白基因的高表达可能表示红细胞污染，这可能不相关于所研究的细胞类型。

dim(sce.all)

## [1] 11548  9580

table(sce.all$orig.ident)

## 
##    1    2    3    4    5    6 
##  350  560  812  671 5995 1192

 3.降维聚类分群
 
   使用Seurat包和harmony包对单细胞RNA测序数据进行标准化、降维、批次校正和聚类分析，下面这段代码意思是：识别可变基因并进行数据标准化；使用PCA对数据进行初步降维；使用Harmony方法进行批次校正；在经过校正后的降维空间中进行聚类分析；使用UMAP和t-SNE两种非线性降维技术进行可视化。

f = "obj1.Rdata"
library(harmony)

if(!file.exists(f)){
  sce.all = sce.all %>% 
  NormalizeData() %>%  #数据标准化
  ScaleData() %>%   #数据缩放
  FindVariableFeatures(selection.method = "vst", nfeatures = 3000) %>%   #寻找可变特征
  RunPCA(pc.genes = [email protected])  %>%  #主成分分析（PCA）
  RunHarmony("orig.ident") %>%  #批次校正
  FindNeighbors(dims = 1:15, reduction = "harmony") %>%  #找到每个细胞的邻居细胞
  FindClusters(resolution = 1) %>%  #聚类分析
  RunUMAP(dims = 1:15,reduction = "harmony") %>%  #使用经过Harmony校正后的前15个主成分进行UMAP降维
  RunTSNE(dims = 1:15,reduction = "harmony")
  save(sce.all,file = f)
}
load(f)

ElbowPlot(sce.all)

     ElbowPlot函数绘制主成分分析（PCA）中的“肘部图”。肘部图用于确定在降维分析中选择多少个主成分是合适的。

     在肘部图中，x轴表示主成分的数量，y轴表示每个主成分的标准差。通常，我们选择曲线开始变得平坦的点作为主成分的数量，即“肘部”位置。这个位置显示了添加更多主成分对解释数据变异的重要性开始减弱的点。

#可视化降维后的数据，通常用于展示聚类结果
p1 <- DimPlot(sce.all, reduction = "umap",label = T)+NoLegend();p1

4.计算marker基因

    下面这段代码用来寻找并保存所有簇（cluster）的标志基因，并从中挑选出每个簇中表达水平最高的两个基因。

f = "markers.Rdata"
if(!file.exists(f)){
  allmarkers <- FindAllMarkers(sce.all,  #寻找所有簇的标志基因
                               only.pos = TRUE, #仅返回正向标志基因（即在簇中上调的基因）
                               min.pct = 0.25, #基因在至少25%的细胞中表达
                               logfc.threshold = 0.25) #对数fold变化阈值为0.25，筛选出表达水平显著上调的基因
  save(allmarkers,file = f)
}
load(f)

mks = allmarkers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
g = mks$gene
g

##  [1] "CCL18"        "SDC3"         "LINC00842"    "TMEM151A"     "LIPA"        
##  [6] "IGSF6"        "PILRA"        "CD86"         "FABP5"        "RPS26"       
## [11] "CXCL3"        "IL8"          "S100A12"      "C19orf59"     "CLDN9"       
## [16] "CXCL17"       "SRPX"         "GRIA3"        "CD3E"         "XCL1"        
## [21] "HS3ST1"       "FABP4"        "RP11-16E12.2" "IGLL1"        "KIF20A"      
## [26] "ANLN"         "LYPD6B"       "FOXJ1"        "BFSP2"        "TNFRSF13B"   
## [31] "PLD4"         "CD1E"         "KIAA0101"     "CDK1"         "CCL19"       
## [36] "BCL2L14"      "SPINK2"       "CDCA7"

5.注释亚群 

    手动注释

#定义标志基因列表
markers = list(epithelial = "EPCAM",
               macrophages = c("CD14","AIF1","CSF1R","CD163"),
               CAFs = c("PDPN","DCN","THY1"),
               DC = c("CD1C","CD1E","CCR7","CD83"),
               Bcells = c("CD19","CD79A","CD79B"),
               Tcells = c("CD2","CD3D","CD3E","CD3G"),
               erythrocytes = "GATA1")

#DotPlot函数用于绘制标志基因在单细胞数据中的表达情况。每个点的大小表示该基因在某个簇中的表达比例，颜色表示表达水平

DotPlot(sce.all,features = markers)+
  RotatedAxis() #将x轴标签旋转

celltype = read.table("celltype.txt") #根据DotPlot图手动填写，以标识每个细胞簇对应的细胞类型
celltype

##    V1          V2
## 1   0 macrophages
## 2   1        CAFs
## 3   2 macrophages
## 4   3 macrophages
## 5   4 macrophages
## 6   5 macrophages
## 7   6 macrophages
## 8   7  epithelial
## 9   8        CAFS
## 10  9           T
## 11 10 macrophages
## 12 11           B
## 13 12        CAFs
## 14 13  epithelial
## 15 14           B
## 16 15 macrophages
## 17 16 macrophages
## 18 17          DC
## 19 18           T

 将注释的细胞类型画到图上去

new.cluster.ids <- celltype$V2
names(new.cluster.ids) <- levels(sce.all)
new.cluster.ids

##             0             1             2             3             4 
## "macrophages"        "CAFs" "macrophages" "macrophages" "macrophages" 
##             5             6             7             8             9 
## "macrophages" "macrophages"  "epithelial"        "CAFS"           "T" 
##            10            11            12            13            14 
## "macrophages"           "B"        "CAFs"  "epithelial"           "B" 
##            15            16            17            18 
## "macrophages" "macrophages"          "DC"           "T"

#RenameIdents函数用于重命名Seurat对象seu.obj中的簇标识符
seu.obj <- RenameIdents(sce.all, new.cluster.ids)
save(seu.obj,file = "seu.obj.Rdata")

#绘制UMAP降维图
p1 <- DimPlot(seu.obj, 
        reduction = "umap", 
        label = TRUE,  #在图中标注每个簇的标签
        pt.size = 0.5) + NoLegend() #移除图例
p1

也可以自己提取数据，用ggplot2画

dat = as.data.frame(seu.obj@[email protected])
dat$seurat_annotation = [email protected]

ggplot(dat,aes(umap_1,umap_2))+
  geom_point(aes(color = seurat_annotation),size = 0.5)+
  scale_color_brewer(palette = "Set1")+
  theme_classic()

6.样本细胞比例条形图

seu.obj$seurat_annotation = [email protected]

ggplot([email protected],aes(patient,
                             fill = seurat_annotation))+
  geom_bar(position = "fill", alpha = 0.9)+
  scale_fill_brewer(palette = "Set1")+
  theme_classic()

7.小提琴图-展示marker基因

    7.1 marker基因

g = unique(g) #标志基因去重

vln.df <- seu.obj@assays$RNA$scale.data %>%
  t() %>%
  as.data.frame()%>%
  select(g) %>% 
  rownames_to_column("CB") %>% 
  mutate(cluster = seu.obj$RNA_snn_res.1)%>%
  pivot_longer(cols = 2:(ncol(.)-1),
               names_to = "gene",
               values_to = "exp") %>% 
  mutate(gene = factor(gene,levels = g))

head(vln.df)

## # A tibble: 6 × 4
##   CB    cluster gene         exp
##   <chr> <fct>   <fct>      <dbl>
## 1 10x_1 0       CCL18      1.75 
## 2 10x_1 0       SDC3       2.16 
## 3 10x_1 0       LINC00842 -0.254
## 4 10x_1 0       TMEM151A  -0.249
## 5 10x_1 0       LIPA      -0.864
## 6 10x_1 0       IGSF6     -0.468

# 自定义颜色
library(paletteer)
my_color = paletteer_d(`"ggsci::default_nejm"`)
my_color = colorRampPalette(my_color)(length(unique(vln.df$cluster)))

# 画图
p1 <- ggplot(vln.df,aes(cluster,exp),color=factor(cluster))+
  geom_violin(aes(fill=cluster),scale = "width")+
  scale_fill_manual(values = my_color)+
  facet_grid(gene~.,scales = "free_y")+
  scale_y_continuous(expand = c(0,0))+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    axis.title.x.bottom = element_blank(),
    axis.ticks.x.bottom = element_blank(),
    axis.text.x.bottom = element_text(hjust = 1,vjust = NULL,color = "black",size = 14),
    axis.title.y.left = element_blank(),
    axis.ticks.y.left = element_blank(),
    axis.text.y.left = element_blank(),
    legend.position = "none",
    panel.spacing.y = unit(0, "cm"),
    strip.text.y = element_text(angle=0,size = 14,hjust = 0),
    strip.background.y = element_blank()
  )

p1

[ixxmu]

小提琴优化