这里是佳奥！单细胞转录组的学习并不轻松，进度过半让我们继续吧。

1 再次了解文章数据

rm(list = ls()) 
Sys.setenv(R_MAX_NUM_DLLS=999)
options(stringsAsFactors = F)
## 首先载入文章的数据
load(file='../input.Rdata')
load(file='../input_rpkm.Rdata')
counts=a
# using raw counts is the easiest way to process data through Seurat.
counts[1:4,1:4];dim(counts)
dat[1:4,1:4];dim(dat)
dat=log2(dat+1)
library(stringr) 
meta=df
head(meta) 
gs=read.table('top18-genes-in-4-subgroup.txt')[,1]
gs
library(pheatmap)
pheatmap(dat[gs,])
pheatmap(dat[gs,],cluster_rows = F)
tmp=data.frame(g=meta$g)
rownames(tmp)=colnames(dat)
pheatmap::pheatmap(dat[gs,],annotation_col = tmp)

table(apply(dat,2,function(x) sum(x>0) )>2000)
FALSE  TRUE 
   76   692 
            
table(apply(dat,1,function(x) sum(x>0) )>4)
 TRUE 
12689 
            
dat=dat[apply(dat,1,function(x) sum(x>0) )>4,
        apply(dat,2,function(x) sum(x>0) )>2000]    

QQ截图20220902142358.png

2 使用Seurat Package（新版参考文档修改）

新版代码参考：

https://cloud.tencent.com/developer/article/1606525

rm(list = ls()) 
Sys.setenv(R_MAX_NUM_DLLS=999)

## 首先载入文章的数据
load(file='../input.Rdata')
counts=a

# using raw counts is the easiest way to process data through Seurat.
counts[1:4,1:4];dim(counts)
library(stringr) 
meta=df
head(meta) 

# 下面的基因是文章作者给出的
gs=read.table('top18-genes-in-4-subgroup.txt')[,1]
gs
library(pheatmap)

fivenum(apply(counts,1,function(x) sum(x>0) ))
boxplot(apply(counts,1,function(x) sum(x>0) ))
fivenum(apply(counts,2,function(x) sum(x>0) ))
hist(apply(counts,2,function(x) sum(x>0) ))
# 上面检测了 counts 和 meta 两个变量，后面需要使用
           
fivenum(apply(counts,1,function(x) sum(x>0) ))
0610010B08Rik        Slc2a5        Nhlrc1          Rfx1    ERCC-00171 
            0             0            20           219           768 

fivenum(apply(counts,2,function(x) sum(x>0) ))
SS2_15_0049_L24  SS2_15_0048_K3 SS2_15_0048_L14 SS2_15_0049_G21  SS2_15_0049_B4 
          198.0          3482.5          4248.0          5234.5          8124.0

library(Seurat)
# https://satijalab.org/seurat/mca.html
# 构建 Seurat 需要的对象
# 其中 min.cells 和 min.genes 两个参数是经验值
sce <- CreateSeuratObject(counts, 
                          meta.data =meta,
                          min.cells = 5, 
                          min.features = 2000, 
                          project = "sce")
# 参考：https://github.com/satijalab/seurat/issues/668

sce
?seurat
table(apply(counts,2,function(x) sum(x>0) )>2000)
table(apply(counts,1,function(x) sum(x>0) )>4)

table(apply(counts,2,function(x) sum(x>0) )>2000)
FALSE  TRUE 
   75   693 
table(apply(counts,1,function(x) sum(x>0) )>4)
FALSE  TRUE 
10140 14442             
            
## 可以看到上面的过滤参数是如何发挥作用的
head(sce@meta.data) 
dim(sce@meta.data)

##查看一下meta.data有哪些信息(后续取列信息只能从这里取)
head(sce@meta.data, 5)
               orig.ident nCount_RNA nFeature_RNA g plate  n_g all
SS2_15_0048_A3        SS2     128784         3067 1  0048 2624 all
SS2_15_0048_A6        SS2     131208         3037 1  0048 2664 all 

## 默认使用细胞名字字符串的切割第一列来对细胞进行分组
# 所以在这里只有一个SS2分组信息, 我们就增加一个 group.by 参数  
            
##旧版代码
#VlnPlot(object = sce, 
        features = c("nGene", "nUMI" ), 
        group.by = 'plate',
        ncol = 2)
#VlnPlot(object = sce, 
        features.plot = c("nGene", "nUMI" ), 
        group.by = 'g',
        nCol = 2)
##新版代码
VlnPlot(sce,
        features = c("nFeature_RNA", "nCount_RNA"),
        group.by = 'plate',##批次效应
        ncol = 2)
VlnPlot(sce,
        features = c("nFeature_RNA", "nCount_RNA"),
        group.by = 'g',##层次聚类分类
        ncol = 2)            

同样的发现，普通的层次聚类得到的4组，很明显是检测到的基因数量的差异造成的

QQ截图20220902152201.png QQ截图20220902152246.png

给sce对象增加一个属性ERCC，供QC使用

##旧版代码
#ercc.genes <- grep(pattern = "^ERCC-", x = rownames(x = sce@raw.data), value = TRUE)
#percent.ercc <- Matrix::colSums(sce@raw.data[ercc.genes, ]) / Matrix::colSums(sce@raw.data)

##新版代码
sce[["percent.ercc"]] <- PercentageFeatureSet(sce, pattern = "^ERCC-")
VlnPlot(object = sce, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.ercc" ), 
        group.by = 'g',
        ncol = 3)

下面这个图第一张和第三张就明显是负相关：也说明了检测的ERCC占比越多，检测到的基因数目就越少

QQ截图20220902154516.png

一些可视化函数

VlnPlot(sce,group.by = 'plate',c("Gapdh","Bmp3","Brca1","Brca2","nFeature_RNA"))

QQ截图20220902160723.png

看两组基因的相关性

##旧版代码
#GenePlot(object = sce,  gene1 = "nUMI", gene2 = "nGene")
#GenePlot(object = sce,  gene1 = "Brca1", gene2 = "Brca2")

##新版代码
plot1 <- FeatureScatter(sce, feature1 = "nFeature_RNA", feature2 = "nCount_RNA")
plot2 <- FeatureScatter(sce, feature1 = "Brca1", feature2 = "Brca2")
CombinePlots(plots = list(plot1, plot2))

QQ截图20220902160901.png

看两个细胞的相关性

##旧版代码
#CellPlot(sce,sce@cell.names[3], 
         sce@cell.names[4],
         do.ident = FALSE)

##新版代码
CellScatter(sce,
            colnames(sce[['RNA']])[3],
            colnames(sce[['RNA']])[4])

QQ截图20220902161159.png

预处理之归一化

sce <- NormalizeData(object = sce, 
                     normalization.method = "LogNormalize", 
                     scale.factor = 10000)

# 当然其中都是默认参数，于是直接写这种也是可以的
sce <- NormalizeData(sce)

# 检查一下归一化后的数据
> sce[['RNA']][1:3,1:3]
3 x 3 sparse Matrix of class "dgCMatrix"
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610005C13Rik              .              .      .        
0610007P14Rik              .              .      0.6256969
0610009B22Rik              .              .      .  
# 结果将0存储为.(为了节省空间，这是松散矩阵sparse Matrix的一个特点)

鉴定差异基因HVGs（highly variable features）

在细胞与细胞间进行比较，选择表达量差别最大的，FindVariableFeatures函数，会计算一个mean-variance结果

# 这里需要理解  dispersion 值
# This function is unchanged from (Macosko et al.), 
# but new methods for variable gene expression identification are coming soon. 
# We suggest that users set these parameters to mark visual outliers on the dispersion plot, but the exact parameter settings may vary based on the data type, heterogeneity in the sample, and normalization strategy. 

##旧版代码
#sce <- FindVariableGenes(object = sce, 
                         mean.function = ExpMean, 
                         dispersion.function = LogVMR )

##新版代码
sce <- FindVariableFeatures(sce, selection.method = "vst")

可见V3、V4.X（新版）相比V2（旧版）做了较大的改动：

• V3计算mean.function和FastLogVMR均采用了加速的FastExpMean、FastLogVMR模式
• V3横坐标范围设定参数改成：mean.cutoff；纵坐标改成：dispersion.cutoff
• 鉴定差异基因的算法包含三种：vst(默认)、mean.var.plot、dispersion
• vst：首先利用loess对 log(variance) 和log(mean) 拟合一条直线，然后利用观测均值和期望方差对基因表达量进行标准化。最后根据保留最大的标准化的表达量计算方差
• mean.var.plot: 首先利用mean.function和 dispersion.function分别计算每个基因的平均表达量和离散情况；然后根据平均表达量将基因们分散到一定数量（默认是20个）的小区间（bin）中，并且计算每个bin中z-score
• dispersion：挑选最高离散值的基因
• V3默认选择2000个差异基因，检查方法也不同（V3用VariableFeatures(sce)检查，V2用sce@var.genes检查）

标准化

这里去除一些技术误差，例如UMI、ERCC，之后就不需要考虑ERCC的影响了，下面的代码中vars.to.regress就是对一些技术因素的排除

关于scale的操作过程：ScaleData这个函数有两个默认参数：do.scale = TRUE和do.center = TRUE，然后需要输入进行scale/center的基因名（默认是取前面FindVariableFeatures的结果）。scale和center这两个默认参数应该不陌生了，center就是对每个基因在不同细胞的表达量都减去均值；scale就是对每个进行center后的值再除以标准差（就是进行了一个z-score的操作）

##旧版代码
#head(sce@meta.data) 
#sce <- ScaleData(object = sce, 
                 vars.to.regress = c("nUMI",'nGene',"percent.ercc" ))

##新版代码
all.genes <- rownames(sce)
length(rownames(sce))
sce <- ScaleData(sce, features = all.genes,
                 vars.to.regress = c("nCount_RNA","percent.ercc" ))
# 结果存储在sce[["RNA"]]@scale.data中
# 后面就不需要考虑ERCC序列了。

实际上，scale函数默认是只针对之前鉴定的HVGs进行标准化（版本3中默认得到2000个HVGs），这样操作的结果中降维和聚类不会被影响，但是只对HVGs基因进行标准化，下面画的热图可能会受到全部基因中某些极值的影响。所以为了热图的结果，还是对所有基因进行归一化比较好。

降维之PCA（后续只保留新版代码）

最常用的PCA方法是一种线性降维，它使用标准化的数据

#默认选择之前鉴定的差异基因（2000个）作为input，但可以使用features进行设置；默认分析的主成分数量是50个；然后输出到屏幕上5个主成分，每个主成分
sce <- RunPCA(sce, features = VariableFeatures(object = sce),
              ndims.print = 1:5, nfeatures.print = 5)

# 这50个主成分的全部结果在 sce@reductions$pca@feature.loadings

对print的主成分结果可视化

#它利用的也就是sce@reductions$pca@feature.loadings结果
VizDimLoadings(sce, dims = 1:2, reduction = "pca")

QQ截图20220902162805.png

按批次检查前两个主成分

DimPlot(sce, reduction = "pca",group.by = 'plate')

可以看到，即使前面没有对批次进行校正，这里也没有批次效应

QQ截图20220902162855.png

再看下层次聚类cluster结果在PCA的可视化

# (DimPlot可以通过降维算法选择，默认首选umap，其次tsne，最后pca，来替代PCAPlot、TSNEPlot、UMAPPlot)
DimPlot(sce, reduction = "pca",group.by = 'g')

# 或者直接 PCAPlot(sce, group.by = 'g')

可以看到几个层次聚类的cluster混在了一起，个人推测可能是：

• PCA默认针对前2000个HVGs进行分析，而层次聚类是针对全部基因做的而分成4群，PCA又是一个线性降维模式，因此分的不是特别开

  
  QQ截图20220902163000.png
  

  然后又探索了使用全部基因做PCA：

sce <- RunPCA(sce, features = all.genes,
              ndims.print = 1:5, nfeatures.print = 5)
PCAPlot(sce, group.by = 'g')

QQ截图20220902163055.png

探索异质性来源—DimHeatmap

每个细胞和基因都根据PCA结果得分进行了排序，默认画前30个基因（nfeatures设置），1个主成分(dims设置)，细胞数没有默认值（cells设置）

DimHeatmap(sce, dims = 1:15, cells = 100, balanced = TRUE)

其中balanced表示正负得分的基因各占一半

QQ截图20220902163149.png
那么降维后选多少个主成分来代表整个数据集？ 最简单使用ElbowPlot(sce)看一看，主要关注”肘部“的PC，它是一个转折点 QQ截图20220902163229.png

细胞聚类

# 先根据ElbowPlot挑选了15个PCs，所以这里dims定义为15个
sce <- FindNeighbors(sce, dims = 1:15)

# 然后使用FindClusters() 进行聚类。其中包含了一个resolution的选项，它会设置一个”间隔“值，该值越大，间隔越大，得到的cluster数量越多
sce <- FindClusters(sce, resolution = 0.4)

## 发现700细胞，设resolution=0.4时，得到Number of communities: 4
> table(sce@meta.data$RNA_snn_res.0.4)
  0   1   2   3 
434 169  52  38 

关于这个resolution参数：分辨率越高，看的越清楚，看到的细节越丰富（cluster越多）；反之，如果分辨率调的很低，结果就看的模模糊糊一大坨

另一种降维方法—非线性降维（tSNE）

sce <- RunTSNE(object = sce, dims.use = 1:15, do.fast = TRUE)
DimPlot(sce,reduction = "tsne",label=T)

QQ截图20220902163420.png

tsne结果

看一下tSNE分类结果和之前层次聚类结果比较：

> table(sce$RNA_snn_res.0.4,sce$g)
   
      1   2   3   4
  0 184 249   0   1
  1  40   2 109  18
  2   4  48   0   0
  3   9   1  12  16
# 左侧0-3是tSNE结果，顶部1-4是hclust结果：hclust的第1组对应tsne的第0组；hclust的第2组还是对应tsne的第0组；hclust的最后一组没有完全分开，说明它们差异还是很大的

找marker基因

# 记住这里tsne的分组是0、1、2、3这四组，因此下面取得ident.1 = 1实际上是第2组
# min.pct的意思是：一个基因在任何两群细胞中的占比最低不能低于多少

cluster2.markers <- FindMarkers(sce, ident.1 = 1, min.pct = 0.25)
> head(cluster2.markers, n = 5)
               p_val avg_log2FC pct.1 pct.2    p_val_adj
Lsp1    2.304267e-89   1.775161 0.828 0.084 3.319528e-85
Pdgfra  3.459284e-86   2.685511 0.870 0.151 4.983444e-82
Svep1   3.087089e-85   2.706946 0.840 0.130 4.447261e-81
Tspan11 3.860428e-83   1.264942 0.757 0.065 5.561333e-79
Mfap5   8.768609e-83   4.011682 0.911 0.242 1.263206e-78

对找到的marker基因进行可视化

• VlnPlot：用小提琴图对某些基因在全部cluster中表达量进行绘制

markers_genes =  rownames(head(cluster2.markers, n = 5)) 
VlnPlot(sce, features = markers_genes)

QQ截图20220902163738.png

第二组的marker基因

• FeaturePlot：最常用的可视化=》将基因表达量投射到降维聚类结果中

FeaturePlot(sce, features = markers_genes)

QQ截图20220902163819.png

第2群marker基因映射结果

可视化文献作者给出的基因

会了基本操作以后，可以将文章中的4个细胞亚群的marker基因拿过来，看看它们分别在我们自己结果中的那一组

取原文的72个marker基因。

# 第一步：鼠标复制(前提是pdf中的图片可以选择文字)，
# 因为基因数太多，这里就不全展示了，只展示一部分
tmp <- c('Atp1b2 Notch3 Ano1 Gm13889 Des Aoc3 Ndu fa4l2 Gucy1a3 Esam Gdpd3 Mcam Higd1b Cpe Kcnj8 Abcc9 Rgs4 Sparcl1 Rgs5 Smoc2 Itgbl1 Fbln1 Cdh11 Crabp1 Pdgf ra Svep1 Pdpn Lsp1 Cpxm1 Lrrc15 Cilp Dcn Lum Mfap5 Fbln2 Olfml3 Rnase4 Mki67 Anln Cdca3 2810417H13Rik Tpx2 Cdca8 Fam64a Nuf2 Birc5 Cep55 Ska1 Kif15 Ttk Melk Top2a Pbk Ccna2 Spc25 Tfap2b Col4a6 Tfap2c Eps8l2 Extl1 Aim1 Irx1 Gata3 Col9a1 Col4a5 Chad Smoc1 Col9a3 Scrg1 Mia Cd24a Perp Trf')

# 第二步：这样粘贴过来的一定存在换行符，就像这样：
> tmp
[1] "Atp1b2 Notch3 Ano1 Gm13889 Des\nAoc3\nNdu fa4l2 Gucy1a3\nEsam\nGdpd3\nMcam\nHigd1b\nCpe\nKcnj8\nAbcc9\nRgs4\nSparcl1\nRgs5\nSmoc2\nItgbl1\nFbln1\nCdh11\nCrabp1\nPdgf ra\nSvep1\nPdpn\nLsp1\nCpxm1\nLrrc15\nCilp\nDcn\nLum\nMfap5\nFbln2\nOlfml3\nRnase4\nMki67\nAnln\nCdca3 2810417H13Rik Tpx2\nCdca8\nFam64a\nNuf2\nBirc5\nCep55\nSka1\nKif15\nTtk\nMelk\nTop2a\nPbk\nCcna2\nSpc25\nTfap2b\nCol4a6\nTfap2c\nEps8l2\nExtl1\nAim1\nIrx1\nGata3\nCol9a1\nCol4a5\nChad\nSmoc1\nCol9a3\nScrg1\nMia\nCd24a\nPerp\nTrf"
# 那么就需要先将换行符换成空格
tmp <- gsub("\n"," ",tmp)
> tmp
[1] "Atp1b2 Notch3 Ano1 Gm13889 Des Aoc3 Ndu fa4l2 Gucy1a3 Esam Gdpd3 Mcam Higd1b Cpe Kcnj8 Abcc9 Rgs4 Sparcl1 Rgs5 Smoc2 Itgbl1 Fbln1 Cdh11 Crabp1 Pdgf ra Svep1 Pdpn Lsp1 Cpxm1 Lrrc15 Cilp Dcn Lum Mfap5 Fbln2 Olfml3 Rnase4 Mki67 Anln Cdca3 2810417H13Rik Tpx2 Cdca8 Fam64a Nuf2 Birc5 Cep55 Ska1 Kif15 Ttk Melk Top2a Pbk Ccna2 Spc25 Tfap2b Col4a6 Tfap2c Eps8l2 Extl1 Aim1 Irx1 Gata3 Col9a1 Col4a5 Chad Smoc1 Col9a3 Scrg1 Mia Cd24a Perp Trf"

# 第三步：分割字符串
library(stringr)
paper_marker <- as.character(str_split(tmp,' ',simplify = T))

# 第四步：检查错误
> length(paper_marker)
[1] 74
# 这个不对，原文只有4组*18=72个基因，多了两个，应该是复制粘贴过来出的错误(因为这里小鼠基因名都是首字母大写，于是先找到基因名首字母不是大写的，再将它们替换掉)

# 其中一个就是‘Ndufa4l2’本来是一个，但是粘贴过来成了两个：'Ndu','fa4l2'
paper_marker_1 <- str_replace(paper_marker,c('Ndu','fa4l2'),"Ndufa4l2") 

# 每次都要检查
> setdiff(paper_marker_1,paper_marker)
[1] "Ndufa4l2"

# 另一个就是"Pdgf" 和"ra"    
paper_marker_2 <- str_replace(paper_marker_1,c('Pdgf','ra'),"Pdgfra")

# 每次修改都要检查：
> setdiff(paper_marker_2,paper_marker_1)
[1] "CPdgfrabp1" "Pdgfra"   
# 发现paper_marker_2中由于替换了一个简单的字符"ra"，所以原来的"Crabp1"也被替换了，修改回来即可
paper_marker_2[paper_marker_2=='CPdgfrabp1'] <- 'Crabp1'

paper_marker <- unique(paper_marker_2)
> length(paper_marker)
[1] 72

使用上面得到的基因进行FeaturePlot

# 原文中的第一群(18个基因中，我们只取前6个基因看一下) 
DimPlot(sce,reduction = "tsne",label=T)
FeaturePlot(sce, features = paper_marker[1:6])

结合我们上面tsne的结果看：

QQ截图20220902164444.png
发现原文中的第一群细胞在我们这里也是大体分到了第一群，当然不是很完美的再现，因为发现这些基因还有少量映射到了旁边的第3群、第4群。

# 同理对其他三群
FeaturePlot(sce, features = paper_marker[19:24])
FeaturePlot(sce, features = paper_marker[37:42])
FeaturePlot(sce, features = paper_marker[55:60])

再回到自己分析的结果，找全部的marker

sce.markers <- FindAllMarkers(sce, only.pos = TRUE, 
                              min.pct = 0.25, logfc.threshold = 0.25)

# 这一步过滤好好理解(进行了分类、排序、挑每个cluster的前20个)
library(dplyr)
top20 <- sce.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)

# 看看和原文挑选的基因重合度(约50%重合度，还不错)
intersect(top20$gene,paper_marker)
 [1] "Rgs5"    "Sparcl1" "Rgs4"    "Atp1b2"  "Abcc9"   "Kcnj8"   "Cpe"     "Des"    
 [9] "Ano1"    "Esam"    "Gucy1a3" "Mfap5"   "Itgbl1"  "Smoc2"   "Cpxm1"   "Dcn"    
[17] "Lum"     "Rnase4"  "Cdh11"   "Fbln1"   "Olfml3"  "Lrrc15"  "Fbln2"   "Crabp1" 
[25] "Cilp"    "Top2a"   "Pbk"     "Mki67"   "Spc25"   "Anln"    "Ccna2"   "Col4a5" 
[33] "Col9a3"  "Trf"    

画自己数据的top20基因热图：

DoHeatmap(sce, features = top20$gene) + NoLegend()

QQ截图20220902165354.png

这一次我们顺利的使用新版代码走完了整个流程。

下一篇我们使用另外两个包。

我们下一篇再见！