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

【参考链接】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()

提取数据

#提取colDatacolData(example_sce)#提取表达量的数据矩阵assay（example_sce,“counts”）#提取降维的矩阵reducedDim(example_sce,"PCA")#提取属性值colnames(colData(example_sce))

存储数据

#增加logcount的数据矩阵example_sce <- logNormCounts(example_sce)#增加PCA的降维矩阵example_sce <- runPCA(example_sce, name="PCA2",ncomponents=10)#增加tSNE的降维矩阵example_sce <- runTSNE(example_sce, perplexity=50,dimred="PCA", n_dimred=10)

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

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

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

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] 22354a<-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) #查看变异矩阵

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

三、查看表达值

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

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

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

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

从下方的这张图中，可以看出基因在不同类型的细胞中的差异表达的特异性。
其中比较直观的可以看出interneurons、pyramidal SS这两种类型的细胞中表达量尤其高。
分类协变量，产生分组小提琴。

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

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

Jam2基因与counts矩阵中前四个基因之间的表达协同性的判断。
图中，每一个圆点代表一个细胞。
从图中可以大概的看到Adamts基因与Jam2基因之间的相关性较低。
连续协变量将在每个面板中生成散点图。

四、降维

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")

使用tSNE降维的过程中，引入了PCA的数据，这里的原理并不是很懂。而且还设置了种子以及perplexity的值。

使用tsne分类，可以看到可以按照组织类型基本分开。
PCA可能还会有一些重合的地方。

注明：

在学习的过程中，发现在使用这个包的过程中，用了很多的内置函数。就是一步处理，这可能会有些傻瓜式。所以在设置参数的时候，要明白其意义。以及最终的生物学意义。

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

Seurat 单细胞转录组测序数据分析教程（二）——python(scanpy)
Seurat 单细胞转录组测序数据分析教程(二)--python(scanpy) 文章参考至scanpy官网,做了一个更详细的解读. 数据由来自健康捐赠者的 3k PBMC组成,可从 10x Geno ...
重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述）
原文链接: https://www.embopress.org/doi/10.15252/msb.20188746 主编评语这篇文章最好的地方不只在于推荐了工具,提供了一套分析流程,更在于详细介绍了 ...
包教会一对一跟着CNS学单细胞测序（含空间转录组、chipseq、RNAseq、Atacseq 和外显子等）3月13日开始...
报名成功立马送您往期所有视频预习本班实行"包教包会.一对一指导服务",即如果报本班,不仅有同步回放视频,而且一对一指导服务,解决学完无法消化问题.学不会免费继续学,直到学会为止. ...
单细胞测序流程（三）质控和数据过滤——Seurat包分析，小提琴图和基因离差散点图
质控和数据过滤准备工具:R. 准备数据:上期经过整理的数据geneMatrix. 注意事项:R的安装目录和文件所在位置都不可有英文. R 语言所需安装的包: #if (!requireNamespa ...
如何在Science、Nature等国际顶刊发文，分子对接、深度学习基因组学，分子动力学、单细胞测序复现文章
基因组学(genomics)是对生物体所有基因进行集体表征.定量研究及不同基因组比较研究的一门交叉生物学学科,基因组学的目的是对一个生物体所有基因进行集体表征和量化,并研究它们之间的相互关系及对生物体 ...
IF：8+ 单细胞测序揭示肝细胞癌的免疫抑制概况
点击关注,桓峰基因桓峰基因的教程不但教您怎么使用,还会定期分析一些相关的文章,学会教程只是基础,但是如果把分析结果整合到文章里面才是目的,觉得我们这些教程还不错,并且您按照我们的教程分析出来不错的结 ...
Nature methods | Alevin-fry，一种高效准确的单细胞测序数据预处理工具
随着单细胞以及单核测序(single-cell and single-nucleus RNA-sequencing)的快速发展以及逐渐普及,越来越多的单细胞测序数据集在近几年不断的出现.这些数据集不仅 ...
哇！单细胞测序-配体受体互作分析原来可以这么简单又高大上！
NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅综述:三万字长文读 ...
单细胞测序分析之小技巧之for循环批量处理数据和出图
"harmony"整合不同平台的单细胞数据之旅生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChI ...
新冠患者样本单细胞测序文献汇总
科研工作者的信仰就是将真相大白于天下 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单 ...

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

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

最新文章

热门文章