【单细胞转录组 实战】三、rawCounts和rpkmNorma
2022-08-31 本文已影响0人
佳奥
这里是佳奥！我们开始解析作者提供的表达矩阵吧！
本次的代码在这里：
https://github.com/jmzeng1314/scRNA_smart_seq2
下载.zip后解压，开始下游分析！
1 安装R包
rm(list = ls())
Sys.setenv(R_MAX_NUM_DLLS=999) ##Sys.setenv修改环境设置，R的namespace是有上限的，如果导入包时超过这个上次就会报错,R_MAX_NUM_DLLS可以修改这个上限
options(stringsAsFactors = F) ##options:允许用户对工作空间进行全局设置，stringsAsFactors防止R自动把字符串string的列辨认成factor
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager") ##判断是否存在BiocManager包，不存在的话安装
library(BiocManager)
BiocManager::install(c( 'scran'),ask = F,update = F)
BiocManager::install("TxDb.Mmusculus.UCSC.mm10.knownGene",ask = F,update = F)
BiocManager::install("org.Mm.eg.db",ask = F,update = F)
BiocManager::install("genefu",ask = F,update = F)
BiocManager::install("org.Hs.eg.db",ask = F,update = F)
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene",ask = F,update = F)
install.packages("ggfortify")
install.packages("FactoMineR")
install.packages("factoextra")
2 解析表达矩阵
上面的是rawCounts表达矩阵，下面的是Normalized归一化后的表达矩阵。
rpkm：单纯比较基因reads数量是没有意义的（测序量不一样，基因长度不一样，测序文库大小不一样），所以需要归一化。
QQ截图20220831102110.png
2.1 rawCounts矩阵
##载入表达矩阵
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',
header = T ,sep = '\t') ##把表达矩阵文件载入R，header=T :保留文件头部信息，seq='\t'以tap为分隔符
##每次都要检测数据
a[1:6,1:4] #对于a矩阵取第1~6行，第1~4列
##读取RNA-seq的 counts 定量结果，表达矩阵需要进行简单的过滤
dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),]
#筛选表达量合格的行(基因), 列（细胞）数不变
##ncol()返回矩阵的列数值；floor()四舍五入取整数；
##function() 定义一个函数；sum()求和
#上面的apply()指令代表对矩阵a进行行计算，判断每行表达量>1的样本总个数，并筛选出细胞表达量合格的基因（行）
#第一个参数是指要参与计算的矩阵——a
#第二个参数是指按行计算还是按列计算，1——表示按行计算，2——按列计算；
#第三个参数是指具体的运算参数,定义一个函数x（即表达量为x）
#对每行中x>1的列（即样本数）求和，即得出的是每行中表达量大于1的样本数，
#然后再筛选出大于floor(ncol(a)/50)的行，这样的行（基因）的细胞表达量才算合格
#因为2%的细胞有表达量，所以对于768个细胞样本，每个行(基因)在细胞中的表达至少要有15.36（约等于15）个样本表达才算合格
## 2 % 的细胞有表达量
##这里是演示归一化如何计算的
dat[1:4,1:4]
sum(dat[,3])
# 0610007P14Rik in SS2_15_0048_A5
log2(18*1000000/sum(dat[,3])+1)
## 18 -- > 6.459884 ## SS2_15_0048_A5
##归一化的一种选择，这里是CPM(count-per-million，每百万碱基中每个转录本的count值)
##CPM只对read count相对总reads数做了数量的均一化，去除文库大小差异。
dat=log2(edgeR::cpm(dat)+1)
dat[1:4,1:4]
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
0610007P14Rik 0 0 6.459884 6.313884
0610009B22Rik 0 0 0.000000 0.000000
0610009L18Rik 0 0 0.000000 0.000000
0610009O20Rik 0 0 2.544699 3.025273
##下面介绍一下 dist 函数
x=1:10
y=2*x
z=rnorm(10)
tmp=data.frame(x,y,z)
dist(tmp)
head(tmp)
dist(t(tmp))
cor(tmp)
dist(t(scale(tmp)))
##可以看到dist函数计算样本直接距离和cor函数计算样本直接相关性，是完全不同的概念。虽然我都没有调它们两个函数的默认的参数。
# 总结：
# - dist函数计算行与行（样本）之间的距离
# - cor函数计算列与列（样本）之间的相关性
# - scale函数默认对每一列（样本）内部归一化
# - 计算dist之前，最好是对每一个样本（列）进行scale一下
##层次聚类，近800细胞。
##原始表达矩阵转置后，细胞在行，所以计算的是细胞与细胞之间的距离。
hc=hclust(dist(t(dat))) ##样本间层次聚类
## statquest 有详细讲解背后的统计学原理。
class(hc)
?plot.hclust
## 查看说明书。
plot(hc,labels = FALSE)
#t:矩阵转置，行转列，列转行
#分类时常常需要估算不同样本之间的相似性(Similarity Measurement)
#这时通常采用的方法就是计算样本间“距离”(Distance)。
#dist函数是R语言计算距离的主要函数。dist函数可以计算行与行两两间的距离。
#所以之前的矩阵里面行是基因，转置后行是样本，因为我们要计算样本与样本之间的距离。
#dist()函数计算变量间距离
#hclust函数用来层次聚类
clus = cutree(hc, 4)##对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus)##转换为因子属性
table(group_list)##统计频数
group_list
1 2 3 4
312 300 121 35
##提取批次信息
colnames(dat) #取列名
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3] #取列名，以'_'号分割，提取第三列。
#str_split()函数可以分割字符串
table(plate)
n_g = apply(a,2,function(x) sum(x>1)) #统计每个样本有表达的有多少行（基因）
# 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。
# 而且大部分单细胞转录组技术很烂，通常超过75%的基因都没办法检测到。
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
##样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】
df$all='all' #添加列，列名为"all"，没事意思，就是后面有需要
metadata=df
save(a,dat,df,file = '../input.Rdata')
##保存a,dat,df这变量到上级目录的input.Rdata
##因为另外一个项目也需要使用这个数据集，所以保存到了上级目录。
##查看一下n_g
hist(n_g)
hist(n_g,breaks = 30)
QQ截图20220831111603.png
2.2 rpkmNormalized矩阵
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',header = T ,sep = '\t')
##把表达矩阵文件载入R，header=T :保留文件头部信息，seq='\t'以tap为分隔符
##每次都要检测数据
a[1:6,1:4] #对于a矩阵取第1~6行，第1~4列
## 读取RNA-seq的 counts 定量结果，表达矩阵需要进行简单的过滤
dat=a[apply(a,1, function(x) sum(x>0) > floor(ncol(a)/50)),] #筛选表达量合格的行,列数不变
dat[1:4,1:4]
#层次聚类
hc=hclust(dist(t(log(dat+0.1)))) ##样本间层次聚类
##如果是基因聚类，可以选择 wgcna 等算法
## statquest
plot(hc,labels = F)
clus = cutree(hc, 4) #对hclust()函数的聚类结果进行剪枝，即选择输出指定类别数的系谱聚类结果。
group_list= as.factor(clus) ##转换为因子属性
table(group_list) ##统计频数
group_list
1 2 3 4
344 189 217 18
QQ截图20220831125306.png
##提取批次信息
colnames(dat) #取列名
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3] #取列名，以'_'号分割，提取第三列。
#str_split()函数可以分割字符串
table(plate)
n_g = apply(a,2,function(x) sum(x>0)) #统计每个样本有表达的有多少行（基因）
# 这里我们定义， reads数量大于1的那些基因为有表达，一般来说单细胞转录组过半数的基因是不会表达的。
# 而且大部分单细胞转录组技术很烂，通常超过75%的基因都没办法检测到。
df=data.frame(g=group_list,plate=plate,n_g=n_g) #新建数据框(细胞的属性信息)
##(样本为行名，列分别为：样本分类信息，样本分组，样本表达的基因数【注意：不是表达量的和，而是种类数或者说个数】)
df$all='all' #添加列，列名为"all"，后面有需要
metadata=df
save(dat,metadata,file = '../input_rpkm.Rdata') #保存dat,df这变量到上级目录的input_rpkm.Rdata
##因为另外一个项目也需要使用这个数据集，所以保存到了上级目录。
3 表达矩阵内部异质性
初步处理了表达矩阵，我们开始正式的分析探索。
文章的图都是基于rpkm值来画的。
3.1 rawCounts矩阵
load(file = '../input.Rdata')
a[1:4,1:4]
head(df) #head()函数显示操作前面的信息，默认前6行
## 载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
##注意:变量a是原始的counts矩阵，变量dat是log2CPM后的表达量矩阵。
group_list=df$g #'$'符，取列，取metadata矩阵的g列,取出层级聚类信息
table(group_list) ##这是全部基因集的聚类分组信息
cg=names(tail(sort(apply(dat,1,sd)),100)) ##取表达量标准差最大的100行的行名
#这个前面演示过一次，dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),] #筛选表达量合格的行,列数不变
#对dat矩阵每行求标准差，排序，取最后的100行，并获得后100行的行名（探针名）
#sd()求标准差，对dat矩阵每一行的counts求标准差
#sort()函数,排序;
#tail()函数，显示操作对象后面的信息，默认后6行，这里设定取后100行
#names()函数，获取或设置对象的名称
library(pheatmap)
##画热图,针对top100的sd的基因集的表达矩阵,没有聚类分组
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F,
filename = 'all_cells_top_100_sd.png')
QQ截图20220831143434.png
##继续针对top100的sd的基因集的表达矩阵，归一化，分组画图
##这部分是示例，归一化不改变数据性质，但是改变数据分布（值域）
x=1:10;plot((x))
scale(x);plot(scale(x))
#scale()函数去中心化和标准化
n=t(scale(t(dat[cg,])))
##对每个探针的表达量进行去中心化和标准化
n[n>2]=2 #矩阵n中归一化后，大于2的项，赋值使之等于2（相当于设置了一个上限）
n[n<-2]= -2 #小于-2的项，赋值使之等于-2（相当于设置了一个下限）
n[1:4,1:4]
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
Kitl 0.9165638 1.2684450 1.4703938 0.6289595
Mmp11 -2.0000000 0.3807272 0.5857658 0.1026709
Atad1 0.4275134 -1.0657421 0.8755671 -1.0657421
Btg2 0.7810682 0.9129418 0.7773478 1.4196704
#pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list) #制作细胞（样本）分组矩阵
rownames(ac)=colnames(n) ##ac的行名（样本名）等于n的列名（样本名）
##判断分组矩阵的行（样本数）和表达矩阵的列（样本数）是否相等
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree1.png')
QQ截图20220831144207.png
##针对top100的sd的基因集的表达矩阵 重新进行 聚类并且分组
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n<-2]=-2
n[1:4,1:4]
##这个聚类分组只是对top100的sd的基因集，从这里开始
hc=hclust(dist(t(n)))
clus = cutree(hc, 4)
group_list=as.factor(clus)
table(group_list) ##这个聚类分组信息是针对top100的sd的基因集的，和针对全部基因集的分组结果不一样
table(group_list,df$g) ## 其中 df$g 是前面步骤针对全部表达矩阵的层次聚类结果
##下面针对本次挑选100个基因的表达矩阵的层次聚类结果进行热图展示
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames = F,show_rownames = F,
annotation_col=ac)
pheatmap(n,show_colnames = F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree_2.png')
dev.off() ##关闭画板
##先对整个基因集聚类分组，再对top100的sd的基因集画图
##和选取top100的sd的基因集，聚类分组再画图，结果聚类分组信息不同
##两次的聚类分组信息不同，画出的图不同
QQ截图20220831145237.png
3.2 PCA分析— rawCounts矩阵
##防止下面操作把数值搞坏的一个备份
dat_back=dat
dat=dat_back ##表达矩阵数据
dat[1:4,1:4]
dat=t(dat)
dat=as.data.frame(dat) ##转换为数据框
dat=cbind(dat,group_list ) ##cbind()合并列（横向追加）;添加分组信息
dat[1:4,1:4]
## 表达矩阵可以随心所欲的取行列，基础知识需要打牢。
dat[1:4,12197:12199]
dat[,ncol(dat)] #ncol()列，返回列长值
table(dat$group_list)
library("FactoMineR")
library("factoextra")
# The variable group_list (index = ) is removed
# before PCA analysis
## 这里的PCA分析，被该R包包装成一个简单的函数，复杂的原理后面讲解。
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$group_list, # color by groups 颜色组
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses 集中成椭圆
legend.title = "Groups")
## 事实上还是有很多基因dropout非常严重。
ggsave('all_cells_PCA.png')
QQ截图20220831152217.png
3.3 rpkmNormalized矩阵
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
# a[1:4,1:4]
dat[1:4,1:4]
head(metadata) #head()函数显示操作前面的信息，默认前6行
##载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
##注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。
group_list=metadata$g #'$'符，取列，取metadata矩阵的g列,取出层级聚类信息
table(group_list) ##这是全部基因集的聚类分组信息
cg=names(tail(sort(apply(dat,1,sd)),100)) ##取表达量标准差最大的100行的行名
# 这个前面演示过一次，dat=a[apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50)),] #筛选表达量合格的行,列数不变
#对dat矩阵每行求标准差，排序，取最后的100行，并获得后100行的行名（探针名）
#sd()求标准差，对dat矩阵每一行的counts求标准差
#sort()函数,排序;
#tail()函数，显示操作对象后面的信息，默认后6行，这里设定取后100行
#names()函数，获取或设置对象的名称
library(pheatmap)
##后面要画热图，需要取对数，不然会被最大值覆盖掉
mat=log2(dat[cg,]+0.01)
##画热图,针对top100的sd的基因集的表达矩阵,没有聚类分组
pheatmap(mat,show_colnames =F,show_rownames = F)
pheatmap(mat,show_colnames =F,show_rownames = F,
filename = 'all_cells_top_100_sd.png')
QQ截图20220831152744.png
上图还是难以解释某些基因在细胞中的表达量是高还是低，难以区分（横向看颜色都差不多）。所以我们并不想知道原始表达量的高低，我们想知道该基因在不同样本表达的高低。
##针对top100的sd的基因集的表达矩阵，归一化，分组画图
x=1:10;plot((x))
scale(x);plot(scale(x))
n=t(scale(t( mat )))
#scale()函数去中心化和标准化
##对每个探针的表达量进行去中心化和标准化
n[n>2]=2 #矩阵n中归一化后，大于2的项，赋值使之等于2（相当于设置了一个上限）
n[n<-2]=-2 #小于-2的项，赋值使之等于-2（相当于设置了一个下限）
n[1:4,1:4]
# pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list) #制作细胞（样本）分组矩阵
rownames(ac)=colnames(n) ##ac的行名（样本名）等于n的列名（样本名）
##判断分组矩阵的行（样本数）和表达矩阵的列（样本数）是否相等
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree1.png')
dev.off() ##关闭画板
QQ截图20220831153042.png
重新进行分组
## 针对top100的sd的基因集的表达矩阵 进行重新 聚类并且分组。
n=t(scale(t( mat )))
n[n>2]=2
n[n<-2]= -2
n[1:4,1:4]
##这个聚类分组只是对top100的sd的基因集
hc=hclust(dist(t(n)))
clus = cutree(hc, 4)
group_list=as.factor(clus)
table(group_list) ##这个聚类分组信息是针对top100的sd的基因集的，和针对全部基因集的分组结果不一样
table(group_list,metadata$g) ## 其中 metadata$g 是前面步骤针对全部表达矩阵的层次聚类结果。
## 下面针对本次挑选100个基因的表达矩阵的层次聚类结果进行热图展示。
ac=data.frame(g=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_cutree_2.png')
dev.off() ##关闭画板
QQ截图20220831154233.png
3.4 PCA分析— rpkmNormalized矩阵
dat_back=dat ##防止下面操作把数值搞坏的一个备份
dat=dat_back ##表达矩阵数据
dat[1:4,1:4]
dat=t(dat)
dat=as.data.frame(dat) ##转换为数据框
dat=cbind(dat,group_list ) ##cbind()合并列（横向追加）;添加分组信息
dat[1:4,1:4]
dat[1:4,12197:12199]
dat[,ncol(dat)] #ncol()列，返回列长值
table(dat$group_list)
library("FactoMineR")
library("factoextra")
# The variable group_list (index = ) is removed
# before PCA analysis
## 这里的PCA分析，被该R包包装成一个简单的函数，复杂的原理后面讲解。
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$group_list, # color by groups 颜色组
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses 集中成椭圆
legend.title = "Groups")
## 事实上还是有很多基因dropout非常严重。
ggsave('all_cells_PCA.png')
QQ截图20220831153403.png
检查两个表达矩阵就到这里。
下一篇我们开始复现第一张图——平均表达量以及变异系数相关散点图。
我们下一篇再见！
