查看原文
其他

用Scater包分析文章数据

刘小泽 单细胞天地 2022-06-06


课程笔记




粉丝:有单细胞线上课程吗?

小编:什么? 我们的单细胞转录组分析线上课程已经上线好久了,你们竟然都不知道吗,每篇推文后面的课程推荐没人看的吗,小编已哭晕在厕所

好了,戏演完了,下面郑重介绍下我们的单细胞线上课程:(详情戳下方链接) 

全网第一个单细胞转录组数据分析实战视频教程(预售第二波通知)

课程链接在:http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53 


这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获!


目录



第三单元第九讲:使用Scater包


首先再次了解文章数据

单细胞转录组学习笔记-4-获取Github代码包以及准备工作

单细胞转录组学习笔记-5-熟悉文献作者提供的两个表达矩阵

载入数据,创建对象

1rm(list = ls()) 
2Sys.setenv(R_MAX_NUM_DLLS=999)
3## 首先载入文章的数据
4load(file='../input.Rdata')
5# 原始表达矩阵
6counts=a
7counts[1:4,1:4];dim(counts) # 2万多个基因,768个细胞
8library(stringr) 
9# 样本信息
10meta=df
11> head(meta) 
12               g plate  n_g all
13SS2_15_0048_A3 1  0048 2624 all
14SS2_15_0048_A6 1  0048 2664 all
15SS2_15_0048_A5 2  0048 3319 all
16SS2_15_0048_A4 3  0048 4447 all
17SS2_15_0048_A1 2  0048 4725 all
18SS2_15_0048_A2 3  0048 5263 all

Scater需要利用SingleCellExperiment这个对象(需要注意的是,官方友情提示,在导入对象之前,最好是将表达量数据存为矩阵

1options(warn=-1# 全局关闭warning信息
2suppressMessages(library(scater))
3## 创建 scater 要求的对象
4sce <- SingleCellExperiment(
5  assays = list(counts = as.matrix(counts)), 
6  colData = meta
7)
8
9> sce
10class: SingleCellExperiment 
11dim: 24582 768 
12metadata(0):
13assays(1): counts
14rownames(24582): 0610005C13Rik 0610007P14Rik ... ERCC-00170
15  ERCC-00171
16rowData names(0):
17colnames(768): SS2_15_0048_A3 SS2_15_0048_A6 ... SS2_15_0049_P22
18  SS2_15_0049_P24
19colData names(4): g plate n_g all
20reducedDimNames(0):
21spikeNames(0):

预处理

如果要计算CPM值,之前一直使用log2(edgeR::cpm(dat)+1)进行计算,这个包自己做了一个函数:calculateCPM()

1exprs(sce) <- log2(calculateCPM(sce ) + 1)
2# 计算结果保存在
3> assayNames(sce)
4[1"counts"    "logcounts"
5# 然后取出来计算结果也很简单
6> counts(sce)[1:3,1:3]
7              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
80610005C13Rik              0              0              0
90610007P14Rik              0              0             18
100610009B22Rik              0              0              0
11> logcounts(sce)[1:3,1:3]
12              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
130610005C13Rik              0              0       0.000000
140610007P14Rik              0              0       6.458664
150610009B22Rik              0              0       0.000000

注意:表达矩阵的标准命名中,exprslogcounts是同义词,它是为了和老版本的scater函数兼容

1> identical(exprs(sce), logcounts(sce))
2[1TRUE

拿到基因名,检测是否存在线粒体基因、ERCC spike-in(利用feature信息对细胞进行质控)

1genes=rownames(rowData(sce))
2> length(genes[grepl('^MT-',genes)])
3[10
4> length(genes[grepl('^ERCC-',genes)])
5[192
6# 这也是标准的ERCC,目前赛默飞公司提供的就设置了92个不同长度和GC含量的细菌RNA序列

利用calculateQCMetrics函数进行QC:

1# 官方给定的写法是:
2example_sce <- calculateQCMetrics(example_sce, 
3    feature_controls = list(ERCC = 1:20, mito = 500:1000),
4    cell_controls = list(empty = 1:5, damaged = 31:40))
5# 这里添加feature_controls信息即可
6sce <- calculateQCMetrics(sce, 
7                          feature_controls = list(ERCC = grep('^ERCC',genes)))

之后过滤:

1# 至少在5个细胞中有表达量的基因可以留下
2keep_feature <- rowSums(exprs(sce) > 0) > 5
3> table(keep_feature)
4keep_feature
5FALSE  TRUE 
610427 14155 
7sce <- sce[keep_feature,]
8# 然后看看细胞中总表达量
9boxplot(sce$total_features_by_counts )
10# 根据下图设置总表达量低于2000的细胞被舍去
11keep.n <- sce$total_features_by_counts > 2000
12filtered_sce <- sce[,keep.n]
13# 最终过滤了一万多基因,几十个细胞
14> dim(filtered_sce)
15[114155   693

对过滤的结果可视化

首先调几个基因看看基因表达量和384孔板这个批次有没有关系:

1plotExpression(sce, rownames(sce)[1:6],
2    x = "plate", exprs_values = "logcounts"

另外简单看看高表达基因的占比:

默认显示前50个基因。图中每一行表示一个基因,每个线条(bar)表示这个基因在不同细胞的表达量(可以想象成把基因表达量的箱线图转了一下)。圆圈是每个基因表达量的中位数,并且在图中经过了排序。

1plotHighestExprs(sce, exprs_values = "counts")

图中可以发现top50中很多ERCC的存在,但是看到第一个ERCC在很多样本中占比超过了15%,也就是说那些样本中有超过15%的reads都”浪费“在了这种不相关的外源序列上,减少了自身的比对量

降维1-PCA分析

默认情况下,runPCA会根据500个变化差异最显著的feature的log-count值进行计算,当然这个数量可以通过ntop参数修改。

1sce <- runPCA(sce)
2plotPCA(sce)
3# SingleCellExperiment对象中包含了reducedDims接口,其中存储了细胞降维后的坐标,可以用reducedDim、reducedDims函数获取,而具体降维的名称用reducedDimNames获取
4> reducedDimNames(sce)
5[1"PCA"

这样的结果需要再加上表型信息,才能看出来是否有批次效应

1plotReducedDim(sce, use_dimred = "PCA"
2                shape_by= "plate"
3                colour_by= "g")
PCA结果

可以看到,四组层次聚类结果在PCA中也能够分开,并且每组细胞都含有两个细胞板信息,细胞板混杂在一起说明没有批次效应

降维2-tSNE

注意:tsne只是一种降维的方法,最后只给出一个坐标。根据坐标是能明白这一群和那一群能分得开,但还不能确定这一群就是单独的一组。为了找到这个依据,就是要进行聚类

1set.seed(1000)
2sce <- runTSNE(sce, perplexity=10)
3plotTSNE(sce, 
4         shape_by= "plate"
5         colour_by= "g")

从图中可以看到,之前用层次聚类定义的4群细胞,现在主要有3群(第1群和第2群混在了一起)

tSNE结果

降维3- DiffusionMap

1sce <- runDiffusionMap(sce)
2plotDiffusionMap(sce,  
3                 shape_by= "plate"
4                 colour_by= "g")
5
6# 此时sce也记录了,做过三次降维处理,但每次都不冲突,存储在平行位置
7> sce@reducedDims
8List of length 3
9names(3): PCA TSNE DiffusionMap
DiffusionMap结果

聚类1-K-means

以t-SNE结果为例

1# 提取tsne降维后的主成分坐标,这个坐标就相当于一个新的矩阵,只不过不记录表达量而记录坐标(供后续聚类使用)
2> head(sce@reducedDims$TSNE,3)
3         [,1]      [,2]
4[1,] 13.99815 -17.82133
5[2,] 14.94079 -14.26358
6[3,] 10.63692 -23.44417
7# Kmeans需要的参数主要是:一个数值矩阵,一个自定义的分群数量
8kmeans(sce@reducedDims$TSNE,centers = 4)
9# 它的结果主要包含:
10# Available components:
11[1"cluster"      "centers"      "totss"        "withinss"    
12[5"tot.withinss" "betweenss"    "size"         "iter"        
13[9"ifault" 
14
15# 我们选择第一个”cluster“,就是记录的分组信息
16> head(kmeans(sce@reducedDims$TSNE,centers = 4)$cluster)
17[14 4 4 3 1 3
18# 最后将分组信息添加到sce的表型信息中
19colData(sce)$tSNE_kmeans <- as.character(kmeans(sce@reducedDims$TSNE,
20                                                centers = 4)$cluster)
21plotTSNE(sce,  colour_by = "tSNE_kmeans")

聚类2-层次聚类hclust:

1# 重点:将tsne的结果当成一个矩阵即可,让其中的坐标分成4群就行
2hc=hclust(dist( sce@reducedDims$TSNE ))
3clus = cutree(hc, 4
4colData(sce)$tSNE_hc <-  as.character(clus)
5plotTSNE(sce,  colour_by = "tSNE_hc")

聚类3-SC3包

另外,有个R包SC3 也提供了聚类的算法,是基于K-means的,用一下试试

1library(SC3) # BiocManager::install('SC3')
2sce <- sc3_estimate_k(sce) # 先预估一下聚类亚群
3sce@metadata$sc3$k_estimation # 预估13个亚群
4rowData(sce)$feature_symbol <- rownames(rowData(sce))
5
6# 接下来正式运行,kn参数表示的预估聚类数
7# 我们这里自定义为4组
8kn <- 4 
9start=Sys.time()
10sce <- sc3(sce, ks = kn, biology = TRUE# 运行会很慢
11end=Sys.time()
12(dur=end-start)
13# 总共需要5分钟左右
14
15# 会将聚类结果放入表型信息(sce@colData)中去,默认叫sc3_cluster,这里人为改个名称 
16sc3_cluster="sc3_4_clusters"
17# 最后进行可视化--比较之前tSNE的Kmeans聚类和SC3的聚类的一致性
18sc3_plot_consensus(sce, k = kn, show_pdata = c("tSNE_kmeans",sc3_cluster))
比较tSNE-Kmeans和SC3-Kmeans

再用table函数看看二者分群的相关性:

1> table(colData(sce)$tSNE_kmeans,colData(sce)$sc3_4_clusters)
2
3      1   2   3   4
4  1  14  42  44   6
5  2 194   0   0   0
6  3  36 147   0   7
7  4 278   0   0   0

同理,也能看看SC3结果与hclust的聚类相关性

1sc3_plot_consensus(sce, k = kn, show_pdata = c("g",sc3_cluster))




破解肿瘤细胞的病人特异性之谜(第一个纯粹的单细胞公共数据库挖掘高分文章)

scRNA-seq揭示胰腺导管腺癌的瘤内异质性和恶性进展

单细胞转录组得到的基因集如何看生存效果呢?(不妨试试看GSVA)

体外培养大脑类器官重现大脑皮层的细胞多样性

总说要挖掘单细胞公共数据,你先用起来可以吗?

Seurat2和Seurat3共存的方法

Monocle2 踩坑教程(2)

Monocle2 踩坑教程(1)

差异分析及KEGG注释简介

单细胞测序助力研究宿主细胞和病原体之间的关系



如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程

生信技能树(爆款入门培训课)全国巡讲约你

全国巡讲第14、15站-兰州和贵阳(生信入门课附赠GEO/TCGA数据挖掘视频课)

全网第二个单细胞视频课程预售




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存