单细胞测序之scater包数据分析教程复现

2024-03-30 19:13:33

【参考链接】http://127.0.0.1:10033/library/scater/doc/overview.

https://bioconductor.org/packages/devel/bioc/vignettes/SingleCellExperiment/inst/doc/intro.html#1_Motivation

【要求】：1、能够完整的复现文档中的每一步的分析步骤。

2、弄清楚为什么要这样处理，得到的结果的生物学意义是什么？

一、加载scRNAseq包中的示例数据

#安装scRNAseq包
BiocManager::install(c('scRNAseq'),ask = F,update = F)
#加载scRNAseq库
library(scRNAseq)
#加载测试数据集（关于小鼠大脑研究的一个数据集）
example_sce <- ZeiselBrainData()

示例数据的文献参考链接：

Zeisel A et al. (2015). Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq. Science 347(6226), 1138-42.

数据集的概况：

example_sce

class: SingleCellExperiment  #单细胞测序特有的数据对象
dim: 20006 3005 #表达矩阵有20006行，3005列。即一共测了3005个细胞，20006个基因的表达量。数据量挺大的。
metadata(0): #对于这个参数，我暂时不理解是什么意思
assays(1): counts  #表达矩阵的名字
rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l  #行的名字（基因名），如Tspan12
rowData names(1): featureType  #这一行也不懂
colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12 1772058148_F03  #细胞的名称
colData names(10): tissue group # ... level1class level2class #为啥是10组我暂时也不明白
reducedDimNames(0):  #降维矩阵
altExpNames(2): ERCC repeat  #2个矩阵 指的是除内源基因以外的多种特征类型的测序数据，如这里的作为内参的ERCC（这里它在后续处理分析时的作用，我仍旧有点不太明白）

这里SingleCellExperiment的这样一个对象的存储形式还是挺有意思的。这样一个数据集怎样提取和存储数据，后续会做一个小小的补充。

数据结构模型

str(example_sce)

Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
  ..@ int_elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  ..@ int_colData        :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  ..@ int_metadata       :List of 3
  ..@ rowRanges          :Formal class 'CompressedGRangesList' [package "GenomicRanges"] with 5 slots
  ..@ colData            :Formal class 'DataFrame' [package "S4Vectors"] with 6 slots
  ..@ assays             :Formal class 'SimpleAssays' [package "SummarizedExperiment"] with 1 slot
  ..@ NAMES              : NULL
  ..@ elementMetadata    :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  ..@ metadata           : list()

提取数据
存储数据

二、对数据矩阵进行进一步的质量检测

对上一步得到的数据矩阵进行下一步的筛选。

移除受损细胞——>怎样判断是否受损呢？
移除低表达量的基因——>如果存在会怎样影响数据分析的结果？

library(scater)
example_sce<-addPerCellQC(example_sce,subsets=list(Mito=grep("mt-",rownames(example_sce))))

对上述addPerlCellQC()函数的参数进行进一步的解释。其中example_sce是我们上一步分析所得到的数据矩阵。subsets是一个包含有一个或者多个vector的list类型数据，是我们比较感兴趣的特殊系列，如ERCC序列或者线粒体基因。

上述代码中的grep()函数的作用类似于正则表达式。Mito=grep("mt-",rownames(example_sce))，也就是说查找example_sce数据矩阵的列名中包含有mt字符的字符串，并将其取出，赋值为Mito。

addPerCellQC，perCellQCMetrics等函数的作用是计算该数据矩阵的一些属性，来识别和删除可能存在问题的单元格。如果我们还指定了特殊的子集（subset），还会额外计算这种子集中的属性，保存在结果中的subset的属性中（见下方结果）。
perCellQCMetrics(
  x,
  subsets = NULL,
  percent.top = integer(0),
  ...,
  threshold=?,        #可以人为设定筛选阈值，影响最后得到的detected的值
  flatten = TRUE,
  assay.type = "counts",
  use.altexps = TRUE,
  percent_top = NULL,
  exprs_values = NULL,
  use_altexps = NULL
)
x
A numeric matrix of counts with cells in columns and features in rows.

Alternatively, a SummarizedExperiment or SingleCellExperiment object containing such a matrix.

subsets
A named list containing one or more vectors (a character vector of feature names, a logical vector, or a numeric vector of indices), used to identify interesting feature subsets such as ERCC spike-in transcripts or mitochondrial genes.

查看输出结果。

colnames(colData(example_sce)

发现与之前相比的话，增加了“sum”，“detected”等之后的一些列。

[1] "tissue"                  "group #"                 "total mRNA mol"         
 [4] "well"                    "sex"                     "age"                    
 [7] "diameter"                "cell_id"                 "level1class"            
[10] "level2class"             "sum"                     "detected"               
[13] "subsets_Mito_sum"        "subsets_Mito_detected"   "subsets_Mito_percent"   
[16] "altexps_ERCC_sum"        "altexps_ERCC_detected"   "altexps_ERCC_percent"   
[19] "altexps_repeat_sum"      "altexps_repeat_detected" "altexps_repeat_percent" 
[22] "total"

在这里对sum,detected等参数进一步解释一下。

sum：为每一个细胞得到的所有基因的表达量之和，即测序总量。
detected：高于阈值的观察数，即每一个细胞中基因表达量高于阈值（前面设置的threshold）的基因的数量。

写代码验证一下猜测。

data_count <- assay(example_sce, "counts") #提取example_sce矩阵中的数据集
sum(data_count[,1]) #计算第一列（第一个细胞的所有基因表达量）的和

【输出结果】：
[1] 22354  #与sum第一个值一致（见下）  

example_sce$sum[1]

【输出结果】：
[1] 22354

a<-sum(data_count[,1]>0) #判断第一列每一个基因的表达量大于0的基因的数量
table(a) #计数

【输出结果】：
a
4871 #输出的数量与detected的第一个值一致，也即第一个细胞样本中表达量大于0的基因的数量（见下）
   1 

example_sce$detected[1]

【输出结果】：
[1] 4871

绘制以sum（测序总量）为横坐标，detected（测序总量中合格的序列）为纵坐标的关于数据质量的基本分布图。

plotColData(example_sce, x = "sum", y="detected", colour_by="tissue")

plotColData(example_sce, x = "sum", y="subsets_Mito_percent", other_fields="tissue") + facet_wrap(~tissue)

这张图的绘制其实也可以反映出数据的质量问题。

横坐标是测序总量，纵坐标是在这些测序的reads中，线粒体DNA(包括其他我们关注的对照指标，如RECC的含量)所占的比例。如果对照指标的比例高，而低表达其他基因，说明测序质量越低的细胞或者空白的细胞（不是我们想要的）。图像中的每一个点，代表的就是一个细胞。

再绘制前50个（函数默认值）特征基因的表达量。

plotHighestExprs(example_sce, exprs_values = "counts")

这一步要运行相对较长的时间。

下一步我们计算每一个基因的变异率，然后统计不同的因素对变异率的贡献情况。对于查看批次效应的影响，药物对基因的影响这一方面上，有很大的作用。

example_sce <- logNormCounts(example_sce) #对counts数据矩阵标准化

vars <- getVarianceExplained(example_sce, variables=c("well", "group #", "level1class", "level2class")) #计算变异率 #原理未知？

head(vars) #查看变异矩阵

单细胞测序之scater包数据分析教程复现 — 这个图像还是很“狂放”的，
可以看出计算得到的这个变异率还是受level，group影响蛮大的。
后经过证实，level1class指的是不同的细胞类型，level2class是对细胞类型的进一步的细分，
group所指的意思暂时未知。

三、查看表达值

绘制特定基因在不同类型的细胞中的表达值。

"x=？"这一步相当于是设置“协变量”的类型，分类协变量将产生如上所示的分组的小提琴，每个特征一个面板。相比之下，连续协变量将在每个面板中生成散点图。

plotExpression(example_sce, rownames(example_sce)[1:6],x = "level1class", colour_by="tissue")

这对于进一步检查从差异表达测试，伪时间分析或其他分析中识别出的基因可能特别有用。

绘制两个基因在细胞中的共表达情况。

plotExpression(example_sce, rownames(example_sce)[1:4],x = rownames(example_sce)[9]，colour_by="tissue")

四、降维

1、主成分分析（PCA）

example_sce <- runPCA(example_sce, name="PCA",subset_row=rownames(example_sce),ncomponents=10)

runPCA（）是scater包中计算PCA的函数。计算得到的结果存放在example_sce的reduceDim中。subset_row参数可以从大矩阵中选择我们感兴趣的小数据集来进行PCA分析，而ncomponents则规定主成分分析的维度。

提取数据集。

PCA<-reducedDim(example_sce, "PCA")
head(PCA)

                    PC1         PC2        PC3      PC4       PC5        PC6        PC7       PC8        PC9    
1772071015_C02 -24.28594 -4.17250313  1.0828863 22.83442 -24.60986  -9.575431  1.7009407 -3.822859 -0.4216202
1772071017_G12 -23.28309 -2.15269093 -5.2062234 21.49579 -23.54741  -5.796298  2.5012823 -6.223696  0.9188150
1772071017_A05 -28.18094 -0.41968292  0.5577381 22.71358 -23.93278  -9.710344 -0.3517405 -1.069919 -6.1292706
1772071014_B06 -28.08441 -0.08877604 -1.3251801 23.98954 -24.62188  -9.395149  0.8155739 -4.078635 -0.2833760
1772067065_H06 -30.00220 -2.47170792  4.5791827 21.24768 -23.41268 -12.112318  2.3176382 -4.986054  4.8463578
1772071017_E02 -28.38489 -4.94868385  5.1796914 24.92079 -25.33863  -8.551646  0.5537437 -1.946099 -2.6903798

绘制主成分分析的图。

plotPCA(example_sce, colour_by="tissue")

2、t-SNE降维方法

set.seed(1000)
example_sce <- runTSNE(example_sce, perplexity=50, 
                       dimred="PCA", n_dimred=10)
head(reducedDim(example_sce, "TSNE"))
plotTSNE(example_sce, colour_by = "tissue")

码农公寓

相关文章