【单细胞转录组 实战】四、复现文章figure——PCA、tSNE分析和基因数量差异

这里是佳奥!

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

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/

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

我们下一篇再见!

©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 220,458评论 6 513
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 94,030评论 3 396
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 166,879评论 0 358
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 59,278评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 68,296评论 6 397
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 52,019评论 1 308
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,633评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,541评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 46,068评论 1 319
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 38,181评论 3 340
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,318评论 1 352
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,991评论 5 347
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,670评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,183评论 0 23
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,302评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,655评论 3 375
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,327评论 2 358

推荐阅读更多精彩内容