这里是佳奥！

熟悉了两个表达矩阵后，我们开始复现文章的图。

1 平均表达量与相关性散点图

统计学概念：

在单细胞测序中起始RNA含量越低，技术重复样本的基因表达相关性也越低，因此生成标准化的基因表达谱之后，非常重要的一步是评估技术噪音（technical variability）
比较常见的一种方法是计算基因表达值的变异系数的平方（CV^2）

标准差与平均数的比值称为变异系数，记为C.V(Coefficient of Variance)
变异系数又称“标准差率”，是衡量资料中各观测值变异程度的另一个统计量
当进行两个或多个资料变异程度的比较时，如果度量单位与平均数相同，可以直接利用标准差来比较

平均绝对误差（Mean Absolute Deviation），又叫平均绝对离差
它是是所有单个观测值与算术平均值的偏差的绝对值的平均

数据分析流程：

rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')

#a[1:4,1:4]
head(metadata)
## 载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
#注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。

dat[1:4,1:4]
exprSet=dat

mean_per_gene <- apply(exprSet, 1, mean, na.rm = TRUE) #对表达矩阵每行求均值
sd_per_gene <- apply(exprSet, 1, sd, na.rm = TRUE) #对表达矩阵每行求标准差
mad_per_gene <-  apply(exprSet, 1, mad, na.rm = TRUE) #对表达矩阵每行求绝对中位差，但在这篇文章没有体现

##同样的apply函数，多次出现，请务必学透它！

##构造一个数据框来存放结果。
cv_per_gene <- data.frame(mean = mean_per_gene,
sd = sd_per_gene,
mad=mad_per_gene,
cv = sd_per_gene/mean_per_gene)
rownames(cv_per_gene) <- rownames(exprSet)
head(cv_per_gene)

#pairs(cv_per_gene)
with(cv_per_gene,plot(log10(mean),log10(cv)))
with(cv_per_gene,plot(log10(mean),log10(cv^2)))

cv_per_gene$log10cv2=log10(cv_per_gene$cv^2)
cv_per_gene$log10mean=log10(cv_per_gene$mean)

library(ggpubr)

##过滤一下，防止离差值带偏拟合曲线
cv_per_gene=cv_per_gene[cv_per_gene$log10mean < 4, ]
cv_per_gene=cv_per_gene[cv_per_gene$log10mean > 0, ]

ggscatter(cv_per_gene, x = 'log10mean', y = 'log10cv2',
          color = "black", shape = 16, size = 1, #Points color, shape and size
          xlab = 'log10(mean)RPKM', ylab = "log10(cv^2)",
          add = "loess", #添加拟合曲线
          add.params = list(color = "red",fill = "lightgray"),
          cor.coeff.args = list(method = "spearman"), 
          label.x = 3, label.sep = "\n",
          dot.size = 2,
          ylim=c(-0.5, 3),
          xlim=c(0,4) 
)

QQ截图20220831162429.png

2 聚类算法简介——PCA和tSNE

目的是检验批次效应。

rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df) 
##载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
# 注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。

group_list=df$g #所有数据的聚类分组信息
plate=df$plate #批次信息
table(plate) 
##这个时候需要从多个维度来探索两个不同的plate的单细胞群体是否有明显的差别。
#plate=group_list

最流行的细胞群体是否有明显的差别，肯定是hclust分群，热图展现，PCA，tSNE 等等。

PCA分析

如果想了解PCA分析原理，需要阅读：

https://mp.weixin.qq.com/s/Kw05PWD2m65TZu2Blhnl4w

首先我们使用简单的 prcomp 函数来了解 PCA分析。

set.seed(123456789)
#set.seed()产生随机数
#用于设定随机数种子，一个特定的种子可以产生一个特定的伪随机序列，这个函数的主要目的，
#是让你的模拟能够可重复出现，因为很多时候我们需要取随机数，但这段代码再跑一次的时候，
#结果就不一样了，如果需要重复出现同样的模拟结果的话，就可以用set.seed()

library(pheatmap)
library(Rtsne)
library(ggfortify)
library(mvtnorm)

## 同样的正态分布随机表达矩阵，是无法区分开来。

ng=500 
nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc) #创建正态分布随机矩阵500行，20列
#dim()检索或设置对象的维度
a2=rnorm(ng*nc);dim(a2)=c(ng,nc) #因为是随机创建，这两个矩阵不一样
a3=cbind(a1,a2)

colnames(a3)=c(paste0('cell_01_',1:nc),
               paste0('cell_02_',1:nc)) #添加列名
 #paste()粘贴，
 rownames(a3)=paste('gene_',1:ng,sep = '') #添加行名
 pheatmap(a3)
 dist(a3)
 a3=t(a3);dim(a3) ## PCA分析 
 pca_dat <- prcomp(a3, scale. = TRUE) #prcomp()主成分分析
 p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
 # df=cbind(as.data.frame(a3),group=c(rep('b1',20),rep('b2',20)))
 # autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour ='group')+theme_bw()
 print(p)
 # 可以看到细胞无法被区分开来

QQ截图20220831173612.png

tSNE

set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0) # Run TSNE
tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2") #添加列名
  # group=c(rep('b1',20),rep('b2',20))
  # tsnes=as.data.frame(tsnes)
  # tsnes$group=group
  # ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point(aes(col=group))
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()

QQ截图20220831173301.png

更多例子

同样的正态分布随机表达矩阵，但是其中部分细胞+3，可以区分开来：

ng=500
nc=20
a1=rnorm(ng*nc);dim(a1)=c(ng,nc)
a2=rnorm(ng*nc)+3;dim(a2)=c(ng,nc) 
a3=cbind(a1,a2)
colnames(a3)=c(paste0('cell_01_',1:nc),paste0('cell_02_',1:nc))
rownames(a3)=paste('gene_',1:ng,sep = '')
pheatmap(a3)
a3=t(a3);dim(a3) 

##PCA分析 
pca_dat <- prcomp(a3, scale. = TRUE)
p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
print(p)
##这个时候细胞被区分开，而且是很明显的一个主成分。

##tSNE分析
set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=10,theta=0.0)
tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2")
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()

QQ截图20220831174240.png QQ截图20220831174251.png QQ截图20220831174258.png

不同的正态分布随机表达矩阵，可以区分：

ng=600
nc=200
mu1  = rnorm(ng, mean = 1)
##rnorm()函数产生一系列的随机数，随机数个数，均值和标准差都可以设定
mu2  = rnorm(ng, mean = 5)
a1=rmvnorm(nc,mu1);dim(a1) #rmvnorm随机生成多元正太分布数
a2=rmvnorm(nc,mu2) ;dim(a2)
    
a3=rbind(a1,a2);dim(a3) #rbind()行进行合并，就是行的叠加
rownames(a3)=c(paste0('cell_01_',1:nc),paste0('cell_02_',1:nc))
colnames(a3)=paste('gene_',1:ng,sep = '')
pheatmap(a3)

##PCA分析 
pca_dat <- prcomp(a3, scale. = TRUE)
p=autoplot(pca_dat) + theme_classic() + ggtitle('PCA plot')
print(p)

##tSNE分析 
set.seed(42)
tsne_out <- Rtsne(a3,pca=FALSE,perplexity=30,theta=0.0)
tsnes=tsne_out$Y
colnames(tsnes) <- c("tSNE1", "tSNE2")
ggplot(tsnes, aes(x = tSNE1, y = tSNE2))+ geom_point()

QQ截图20220831174701.png QQ截图20220831174710.png QQ截图20220831174718.png

然后我们可以使用高级R包做真实的分析。

3 rpkm矩阵PCA分析

##下面是画PCA的必须操作，需要看不同做PCA的包的说明书。
dat_back=dat

dat=dat_back
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,plate ) #cbind根据列进行合并，即叠加所有列 #矩阵添加批次信息
dat[1:4,1:4]
table(dat$plate)

## 'princomp' can only be used with more units than variables
## Principal component analysis is underspecified if you have fewer samples than data point. 

# pca_dat =  princomp(t(dat[,-ncol(dat)]))$scores[,1:2]
# pca_dat =  prcomp(t(dat[,-ncol(dat)])) 
 
# plot(pca_dat$rotation[,1:2], t='n')
# colors = rainbow(length(unique(dat$plate)))
# names(colors) = unique(dat$plate)
# text(pca_dat$rotation[ , 1:2], labels=dat$plate,col=colors[dat$plate])

rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')

dat[1:4,1:4]
head(metadata) 
plate=metadata$plate

##下面是画PCA的必须操作，需要看不同做PCA的包的说明书。
dat_back=dat

dat=dat_back
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,plate) #cbind根据列进行合并，即叠加所有列 #矩阵添加批次信息
dat[1:4,1:4]
table(dat$plate)
library("FactoMineR")
library("factoextra") 
# The variable plate (index = ) is removed
# before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)
fviz_pca_ind(dat.pca,#repel =T,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = dat$plate, # color by groups
             #palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
) 
ggsave('all_cells_PCA_by_plate.png') 
## 保存你的画布的图到本地图片文件。

QQ截图20220831175352.png

4 每个细胞检测到的基因数量差异

rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df) 

##载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
#注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。
group_list=df$g
plate=df$plate
table(plate)
# n_g = apply(a,2,function(x) sum(x>1)) #统计每个样本有表达的有多少行（基因）

## 绘图必须强推 ggpubr 

rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
# n_g = apply(a,2,function(x) sum(x>0)) #统计每个样本有表达的有多少行（基因）

dat[1:4,1:4]
df=metadata
head(df) 
library(ggpubr)
ggviolin(df, x = "all", y = "n_g", fill = "all", 
         add = "boxplot", add.params = list(fill = "white")) 

library(ggpubr)
ggviolin(df, x = "plate", y = "n_g", fill = "plate",
         #palette = c("#00AFBB", "#E7B800", "#FC4E07"),
         add = "boxplot", add.params = list(fill = "white")) 

library(ggpubr)
ggviolin(df, x = "g", y = "n_g", fill = "g", 
         add = "boxplot", add.params = list(fill = "white"))  + stat_compare_means()

查看基因分布：n_g(number of gene)

QQ截图20220831180650.png

不同批次的情况：分布基本相同。

QQ截图20220831180702.png

查看不同分组的情况：有明显区别。

QQ截图20220831180714.png

说明这个分群是技术差异，没有检测到，没有生物学意义。这样分群是不对的。

因此文章采用的是DBSCAN。

关于画图优化的网站：

http://sthda.com/english/

下一篇我们补充一些乳腺癌领域的背景知识。

我们下一篇再见！