条条道路通罗马—单细胞分群分析 / 开普饭

课程笔记

粉丝：有单细胞线上课程吗？

小编：什么

？我们的单细胞转录组分析线上课程已经上线好久了，你们竟然都不知道吗，每篇推文后面的课程推荐没人看的吗，小编已哭晕在厕所

好了，戏演完了，下面郑重介绍下我们的单细胞线上课程：（详情戳下方链接）

全网第二个单细胞视频课程预售

这个课程笔记栏目记录了学员们学习单细胞转录组课程的学习笔记

希望大家能有所收获！

作者 | 单细胞天地小编刘小泽

课程链接在：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")

“

marker基因的表达量可视化

”

条条道路通罗马—单细胞分群分析