其他
条条道路通罗马—单细胞分群分析
课程笔记
粉丝：有单细胞线上课程吗？
小编：什么
好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接）
这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记
希望大家能有所收获！
作者 | 单细胞天地小编 刘小泽
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二单元最后一讲第11讲：单细胞转录组的其他R包介绍
这一次的问题是：分析单细胞转录组一定要用R包吗？
之前在 差异分析及可视化 中使用monocle的plot_cell_clusters函数画出了PBMC的第4和第10群两种不同T细胞的差异。那么这个分析一定要用包装好的R包吗？不是的，即使不使用别人做的R包，自己也能利用作图原理画出来
问题的重点不在使用什么包，而是做出这个图的原理是什么
就上面这个图而言，为了分群，首先要有一个表达矩阵，还要设置分群信息（例如上面我们就明确知道要分成两群），然后还要进行PCA的线性降维和tSNE进一步非线性降维。有了这些认识，那么首先使用普通的函数来尝试一下
第一种方式：不使用R包，使用常规函数
先加载表达矩阵和分群信息
1options(warn=-1)
2load('patient1.PBMC_RespD376_for_DEG.Rdata')
3# 其中保存了表达矩阵（count_matrix）和分群的信息（cluster）
4> count_matrix[1:4,1:4]
5 AAACCTGCAACGATGG.3 AAACCTGGTCTCCATC.3 AAACGGGAGCTCCTTC.3 AAAGCAATCATATCGG.3
6RP11-34P13.7 0 0 0 0
7FO538757.2 0 0 0 0
8AP006222.2 0 0 0 0
9RP4-669L17.10 0 0 0 0
10> table(cluster)
11cluster
12 4 10
13636 170
根据分群因子，赋予不同颜色
1# 根据分群因子，赋予不同颜色
2color <- cluster
3levels(color) <- rainbow(2)
4> table(color)
5color
6#FF0000FF #00FFFFFF
7 636 170
对表达矩阵过滤
1# 首先是对基因过滤，根据标准差将那些在细胞中的表达量没有变化的基因去掉
2choosed_count <- count_matrix
3choosed_count <- choosed_count[apply(choosed_count, 1, sd)>0,]
4# 看到过滤了5000多个基因
1# 然后选择变化差异最明显的前1000个基因
2choosed_count <- choosed_count[names(head(sort(apply(choosed_count, 1, sd),decreasing = T),1000)),]
然后进行PCA分析
对行进行操作，因此需要转置，另外进行标准化
1pca_out <- prcomp(t(choosed_count),scale. = T)
2library(ggfortify)
3autoplot(pca_out, col=color) +theme_classic()+ggtitle('PCA plot')
PCA的全部结果可以通过str(pca_out)了解，其中坐标的信息在pca_out$x
1> pca_out$x[1:3,1:3]
2 PC1 PC2 PC3
3AAACCTGCAACGATGG.3 -1.4940046 6.707691 -0.7655250
4AAACCTGGTCTCCATC.3 0.1403445 3.239395 0.2672606
5AAACGGGAGCTCCTTC.3 -1.9968613 3.924992 -0.4669826
再用tSNE继续降维
1library(Rtsne)
2# 利用PCA的前5个主成分
3tsne_out <- Rtsne(pca_out$x[,1:6], perplexity = 10,
4 pca = F, max_iter = 2000,
5 verbose = T)
6# tSNE的坐标结果保存在tsne_out$Y中
7tsnes_cord <- tsne_out$Y
8# ggplot作图
9colnames(tsnes_cord) <- c('tSNE1','tSNE2')
10ggplot(tsnes_cord, aes(x=tSNE1, y = tSNE2)) + geom_point(col=color) + theme_classic()+ggtitle('tSNE plot')
第二种方式：使用Seurat
Seurat V2
1PBMC <- CreateSeuratObject(count_matrix,
2 min.cells = 1, min.features = 0,
3 project = '10x_PBMC')
4PBMC <- AddMetaData(object = PBMC, metadata = apply(count_matrix, 2, sum), col.name = 'nUMI_raw')
5PBMC <- AddMetaData(object = PBMC, metadata = cluster, col.name = 'cluster')
6VlnPlot(PBMC, features.plot = c("nUMI_raw", "nGene"),
7 group.by = 'cluster',nCol = 2)
看到第10群细胞的基因数和UMI数比第4群多很多，也就是说它们本身的“起跑线”就有差异，那么最后的PCA结果差异难道是由于基因和UMI数量导致的？
那么如果我们接下来进行校正，会发生什么？
1PBMC <- ScaleData(PBMC, vars.to.regress = c("nUMI_raw", "nGene"),
2 model.use = 'linear',
3 use.umi = F)
4# 这个FindVariableGenes其实和前面的找top1000 sd基因目的一样，就是为了增加分析的有效性，认为高变化的基因存储的信息更值得关注
5PBMC <- FindVariableGenes(object = PBMC, mean.function = ExpMean,
6 dispersion.function = LogVMR,
7 x.low.cutoff = 0.0125, x.high.cutoff = 3,
8 y.cutoff = 0.5)
9length(PBMC@var.genes)
10PBMC <- RunPCA(object = PBMC, pc.genes = PBMC@var.genes)
11PBMC <- RunTSNE(object = PBMC, dims.use = 1:20)
12TSNEPlot(PBMC, group.by='cluster')
这次结果发现，这两群就分不开了。从这里也反映出一些问题：本文的这个热图真的是由于生物学因素导致的吗？
猜想：可能这两群细胞本身表达的基因数量就不同，就是有一些基因在这群细胞表达，在那一群不表达
Seurat V3
1PBMC_V3 <- CreateSeuratObject(counts = count_matrix,
2 min.cells = 1,
3 min.features = 0,
4 project = "10x_PBMC_V3")
5PBMC_V3 <- AddMetaData(object = PBMC_V3, metadata = apply(count_matrix, 2, sum), col.name = 'nUMI_raw')
6PBMC_V3 <- AddMetaData(object = PBMC_V3, metadata = cluster, col.name = 'cluster')
7VlnPlot(PBMC_V3, features = c("nUMI_raw", "nFeature_RNA"), ncol = 2,group.by = 'cluster')
8# 同样也进行校正一下
9PBMC_V3 <- ScaleData(object = PBMC_V3, vars.to.regress = c("nUMI_raw", "nFeature_RNA"),
10 model.use = 'linear', use.umi = FALSE)
11
12PBMC_V3 <- FindVariableFeatures(object = PBMC_V3, mean.function = ExpMean,
13 dispersion.function = LogVMR,
14 mean.cutoff = c(0.0125,3),
15 dispersion.cutoff = c(0.5,Inf))
16length(VariableFeatures(object = PBMC_V3))
17PBMC_V3 <- RunPCA(PBMC_V3, features = VariableFeatures(object = PBMC_V3))
18pbmc_tsne <- RunTSNE(PBMC_V3, dims = 1:20)
19DimPlot(pbmc_tsne, reduction = "tsne",group.by='cluster')
第三种方式：使用monocle
1library(monocle)
2# 1.表达矩阵
3expr_matrix <- as.matrix(count_matrix)
4# 2.细胞信息
5sample_ann <- data.frame(cells=names(count_matrix),
6 cellType=cluster)
7rownames(sample_ann)<- names(count_matrix)
8# 3.基因信息
9gene_ann <- as.data.frame(rownames(count_matrix))
10rownames(gene_ann)<- rownames(count_matrix)
11colnames(gene_ann)<- "genes"
12# 然后转换为AnnotatedDataFrame对象
13pd <- new("AnnotatedDataFrame",
14 data=sample_ann)
15fd <- new("AnnotatedDataFrame",
16 data=gene_ann)
17# 最后构建CDS对象
18sc_cds <- newCellDataSet(
19 expr_matrix,
20 phenoData = pd,
21 featureData =fd,
22 expressionFamily = negbinomial.size(),
23 lowerDetectionLimit=0.5)
24
25cds=sc_cds
26cds <- detectGenes(cds, min_expr = 1)
27expressed_genes <- row.names(subset(cds@featureData@data,
28 num_cells_expressed >= 1))
29length(expressed_genes)
30cds <- cds[expressed_genes,]
31
32cds <- estimateSizeFactors(cds)
33cds <- estimateDispersions(cds)
34disp_table <- dispersionTable(cds) # 挑有差异的
35unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1) # 挑表达量不太低的
36cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id) # 准备聚类基因名单
37# 进行降维
38cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
39 reduction_method = 'tSNE', verbose = T)
40# 进行聚类
41cds <- clusterCells(cds, num_clusters = 4)
42# Distance cutoff calculated to 1.812595
43plot_cell_clusters(cds, 1, 2, color = "cellType")
第四种方式：使用scater
1suppressMessages(library(scater))
2## 创建 scater 要求的对象
3# 在导入对象之前，最好是将表达量数据存为矩阵
4sce <- SingleCellExperiment(
5 assays = list(counts = as.matrix(count_matrix)),
6 colData = data.frame(cluster)
7)
8# 预处理
9exprs(sce) <- log2(calculateCPM(sce ) + 1)
10# 降维
11sce <- runPCA(sce)
12plotReducedDim(sce, use_dimred = "PCA",
13 colour_by= "cluster")
1set.seed(1000)
2sce <- runTSNE(sce, perplexity=10)
3plotTSNE(sce,
4 colour_by= "cluster")
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
（福州、上海见！）全国巡讲第19-20站（生信入门课加量不加价）
”