scRNA-seq数据分析学习(二)

学习R包scater对scRNA-seq数据进行分析

scRNA-seq数据分析(二)

1
2
rm(list = ls())#clear the environment
options(warn=-1)#turn off warning message globally

认识scater

创建scater要求的对象

1
2
3
4
5
6
7
8
9
10
11
12
library(scRNAseq)
data(fluidigm)
assay(fluidigm)<-assays(fluidigm)$rsem_counts
ct<-floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
table(rownames(ct)==0)
sample_ann<-as.data.frame(colData(fluidigm))#数据集的临床信息
#这里需要将表达矩阵做成我们的scater要求的对象
sce<-SingleCellExperiment(
assays=list(counts=ct),
colData=sample_ann
)

过滤

基因层面的过滤

使用calculateQCMetrics函数作用于sce这个单细胞数据对象后,就可以用rowData(object)查看各个基因各项统计指标:

  • mean_counts:平均表达量
  • log10_mean_counts:归一化 log10-scale
  • pct_dropout_by_counts:该基因丢失率
  • n_cell_by_counts:多少个细胞表达了该基因

上面的指标可以用来过滤,也可以自己计算这些统计学指标
主要是过滤掉低表达量的基因,还有线粒体基因ERCC spike-ins的控制

1
2
3
4
5
6
7
8
9
exprs(sce)<-log2(calculateCPM(sce)+1)
genes<-rownames(rowData(sce))#rowData(object)基因相关统计情况
genes[grepl('^MT-',genes)]
genes[grepl('^ERCC-',genes)]
#比较不幸,这个测试数据里面没有线粒体基因,也没有ERCC序列
sce<-calculateQCMetrics(sce,
feature_controls = list(ERCC=grep('^ERCC',genes)))
#后面的分析都是基于sce这个变量,这个是一个SingleCellExperiment对象
#sce是一个list(S4),但这个list可以当成data.frame来使用,非常方便
1
2
3
4
#查看信息
tmp<-as.data.frame(rowData(sce))
colnames(tmp)
head(tmp)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#目前只过滤掉那些在所有细胞都没有表达的基因
#这个过滤条件可以自行调整,可以看到基因数量大幅减少
keep_feature<-rowSums(exprs(sce)>0)>0
table(keep_feature)
sce<-sce[keep_feature,]
sce
```

<img src="https://blog-image-host.oss-cn-shanghai.aliyuncs.com/gyqblog/scemid.JPG"/>


**细胞层面的过滤**
用`colData(object)`可以查看各个样本统计情况
```R
tmp<-as.data.frame(colData(sce))
colnames(tmp)
#尝试查看一些信息
tf<-sce$total_features_by_counts
boxplot(tf)
fivenum(tf)
#还是那句话,过滤没有统一的标准,视具体情况而定,下面举一个例子进行简单的过滤尝试
tmp$pct_counts_in_top_100_features_endogenous<50
table(tmp$pct_counts_in_top_100_features_endogenous<50)
sce<-sce[,tmp$pct_counts_in_top_100_features_endogenous<50]
sce

这个sce真的是非常方便了,虽然是个list,但是却可以像dataframe一样对其进行处理,
至此,我们过滤掉了9000多个基因,4个样本

数据可视化

1
2
3
4
#展示一些基因在不同细胞分类的表达,这里只做一下演示,一般不这么用,不过筛出来的marker基因可以用这种方式展现其在各组间的差异
plotExpression(sce,rownames(sce)[1:6],
x="Biological_Condition",
exprs_values = "logcounts")

PCA

1
2
3
4
5
6
#还可以可视化细胞距离分布
sce<-runPCA(sce)#这里没有进行任何基因的挑选,就直接进行PCA了,与seurat包不一样
reducedDimNames(sce)
#PCA分布图上添加临床信息,同样发现不同GW细胞之间不能很好的分群
#plotPCA(sce)该数据最原始的PCA图
plotReducedDim(sce,use_dimred = "PCA",colour_by = "Biological_Condition")
1
2
#PCA分布图上面添加表达量信息,一般标marker gene,次方法不常用,这里就以一个基因为例示范一下
plotReducedDim(sce,use_dimred = "PCA",colour_by = rownames(sce)[1],size_by = rownames(sce)[1])
1
2
3
##仅仅是选取前20个PC,分群时细胞分布并没有太大区别
sce<-runPCA(sce,ncomponents = 20)
plotPCA(sce,colour_by="Biological_Condition")
1
2
##可以挑选指定的PC来可视化,这里选5个PC
plotPCA(sce,ncomponents=5,colour_by="Biological_Condition")

tSNE

1
2
3
4
#tSNE可视化
set.seed(1234)
sce<-runTSNE(sce,perplexity = 30)
plotTSNE(sce,colour_by="Biological_Condition")
1
2
3
4
5
6
7
#perplexity的设定值不一样,散点的分布也会有所不同,但趋势都是一致的
#perplexity设置为10
sce<-runTSNE(sce,perplexity = 10,use_dimred = "PCA",n_dimred = 10)
plotTSNE(sce,colour_by="Biological_Condition")
#perplexity设置为20
sce<-runTSNE(sce,perplexity = 20,use_dimred = "PCA",n_dimred = 10)
plotTSNE(sce,colour_by="Biological_Condition")

认识SC3

R包SC3处理scRNAseq内置的数据

1
2
3
4
5
6
7
8
9
10
11
#R包SC3
library(SC3)
sce<-sce_back
sce<-sc3_estimate_k(sce)
metadata(sce)$sc3$k_estimation
rowData(sce)$feature_symbol<-rownames(rowData(sce))
#一步运行sc3的所有分析,相当耗时,数据量大的话就真的相当耗时
#这里kn表示的预估聚类数,数据集是已知的,我们这里就设定为4组,具体数据要具体考虑
kn<-4
sc3_cluster<-"sc3_4_clusters"
sce<-sc3(sce,ks=kn,biology = TRUE)

可视化展示

kn就是聚类数

1
2
#热图:比较先验分类和SC3的聚类的一致性
sc3_plot_consensus(sce,k=kn,show_pdata = c("Biological_Condition",sc3_cluster))
1
2
#热图:展示表达量信息
sc3_plot_expression(sce,k=kn,show_pdata = c("Biological_Condition",sc3_cluster))

SC3包找marker基因

1
2
#热图,展示可能的标记基因
sc3_plot_markers(sce,k=kn,show_pdata =c("Biological_Condition",sc3_cluster) )