这里是佳奥！我们开始解析作者提供的表达矩阵吧！

本次的代码在这里：

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

检查两个表达矩阵就到这里。

下一篇我们开始复现第一张图——平均表达量以及变异系数相关散点图。

我们下一篇再见！