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

2021-02-27 19:02:02 阅读：570 来源： 互联网

标签：sce 测序 sum 基因 detected 复现 example 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()

提取数据
存储数据

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

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

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

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

这个图像还是很“狂放”的，
可以看出计算得到的这个变异率还是受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可能还会有一些重合的地方。

标签：sce,测序,sum,基因,detected,复现,example,scater
来源： https://blog.csdn.net/weixin_40640700/article/details/113955890

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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