单细胞测序之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") 

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

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

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

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

横坐标是测序总量,纵坐标是在这些测序的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")

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

单细胞测序之scater包数据分析教程复现
从下方的这张图中,可以看出基因在不同类型的细胞中的差异表达的特异性。
其中比较直观的可以看出interneurons、pyramidal SS这两种类型的细胞中表达量尤其高。
分类协变量,产生分组小提琴。

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

plotExpression(example_sce, rownames(example_sce)[1:4],x = rownames(example_sce)[9],colour_by="tissue")
单细胞测序之scater包数据分析教程复现
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")

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

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

 

单细胞测序之scater包数据分析教程复现
使用tsne分类,可以看到可以按照组织类型基本分开。
PCA可能还会有一些重合的地方。

 

 

 

上一篇:【hbase】检察HBase集群的一致性并修复


下一篇:深入理解Flask_(美)Jack StoufferPDF高清文档下载