把单细胞表达量矩阵换一个单位

ixxmu commented 3 years ago

https://mp.weixin.qq.com/s/-jAaHAvraGnVLsPLGKTvQA

github-actions[bot] commented 3 years ago

把单细胞表达量矩阵换一个单位 by 生信技能树

一般来说单细胞表达量矩阵都是以基因为单位，我们可以很容易走常规的降维聚类分群并且合理的进行生物学命名，比如我们对官方 pbmc3k 例子，跑标准代码：

library(Seurat)
# https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
## Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")

## Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", 
                           min.cells = 3, min.features = 200)
## Identification of mithocondrial genes
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

## Filtering cells following standard QC criteria.
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & 
                 percent.mt < 5)

## Normalizing the data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", 
                      scale.factor = 10000)

pbmc <- NormalizeData(pbmc)

## Identify the 2000 most highly variable genes
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

## In addition we scale the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), 
               verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")
table(pbmc$seurat_clusters)
# pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25,  logfc.threshold = 0.25, verbose = FALSE)
DimPlot(pbmc, reduction = "umap", group.by = 'seurat_clusters',
        label = TRUE, pt.size = 0.5) 
DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", 
                               "CD14", "FCER1A", "FCGR3A", 
                               "LYZ", "PPBP", "CD8A"),
        group.by = 'seurat_clusters')
## Assigning cell type identity to clusters
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T",
                     "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
 
pbmc$cluster_by_counts=Idents(pbmc)
table(pbmc$cluster_by_counts)

得到的结果如下所示：

降维聚类分群

根据我们的生物学背景，可以进行如下所示的命名：

生物学命名

这个时候的生物学背景因人而异哦，比如我这里对myeloid了解不多，所以就没有细分出来巨噬细胞，单核细胞，树突细胞等等。仅仅是举例而已。

首先可以替换成为转录因子活性

这里推荐一个 R包，DoRothEA on http://bioconductor.org/ and https://github.com/saezlab/dorothea.

获取包自带数据库

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("dorothea")
library(tidyverse)
## We read Dorothea Regulons for Human:
dorothea_regulon_human <- get(data("dorothea_hs", package = "dorothea"))

## We obtain the regulons based on interactions with confidence level A, B and C
regulon <- dorothea_regulon_human %>%
    dplyr::filter(confidence %in% c("A","B","C"))

使用数据库去分析pbmc

## We compute Viper Scores 
library(dorothea)
pbmc <- run_viper(pbmc, regulon,
                  options = list(method = "scale", minsize = 4, 
                                 eset.filter = FALSE, cores = 1, 
                                 verbose = FALSE))

重新进行降维聚类分群：

## We compute the Nearest Neighbours to perform cluster
DefaultAssay(object = pbmc) <- "dorothea"
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc, features = rownames(pbmc), verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)

pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")
 
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()


library(gplots)
balloonplot(table(pbmc$dorothea_snn_res.0.5,
                  pbmc$cluster_by_counts))

同一个10x单细胞数据的两个矩阵各自的降维聚类分群是类似的，但是又不能完全一致：

这个时候的表达量矩阵并不是原始的每个基因在每个细胞里面的counts值了哦，而是每个转录因子在每个细胞的活性情况，我们也可以人工命名，

后面也可以对这样的矩阵（每个转录因子在每个细胞的活性情况），进行差异分析等等，参考；https://bioconductor.org/packages/release/data/experiment/vignettes/dorothea/inst/doc/single_cell_vignette.html

也可以可以替换成为gsea打分

单位就是各个基因集，比如MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据
C7: immunologic signatures: 免疫相关基因集合。

全新的矩阵是基因集在每个细胞的打分矩阵，然后可以再次聚类看看！

作为学徒作业

大家需要拿kegg的约300个通路，对这个pbmc3k单细胞表达量矩阵进行ssGSEA打分后得到矩阵，进行降维聚类分群！

写在文末

没有氛围，一个人吭哧吭哧的完成作业或多或少会遇到困难，如果有热火朝天的讨论，应该是会极大的促进咱们的积极性。所以有必要组建一个群聊，专门的针对于这个学徒作业100个练习题~

老规矩，群聊组建费用18.8元，一个简单的门槛隔绝那些不怀好意的广告营销号！ 前200名可以直接扫描（仍然需要18.8）群聊二维码进群，满200人后我们会统一收款！（每个人都是18.8 元，如果你不同意这个象征性收费，请不要进群哈！）

如果上面的二维码无法进群，说明已经满员了，需要我们生信技能树的官方拉群小助手帮忙拉群哦！！！（名额有限，先到先得！！！）这个时候请直接付款28元给小助手，就可以进群，或者你转发此推文到朋友圈然后截图给小助手，就可以仍然以18.8元进群！

一个简单的门槛，隔绝那些营销号！我们也会在群里共享学徒作业的资料，仅此而已，考虑清楚哦！

长按识别二维码

烦请备注姓名学校单位信息

会在群里共享收集整理好的全部文献

迈向数字游民之路

最近在朋友圈刷到了【一刻talks先见未来大会2020】的演讲视频，其中一个讲者很有意思、就是海尔集团海创汇合伙人，首席生态官——檀林，他的演讲：数字游民时代来临，你做好准备了吗？ 详细演讲视频介绍见：檀林：数字游民时代来临，你做好准备了吗？

让我接触到了“数字游民”这个词，其实最先提出“数字游民”构想的是一位名叫Tsugio Makimoto的日本高管。早在1997年，他便出版了一本名叫《Digital Nomad》的书：“未来的人类社会，高速的无线网络和强大的移动设备会打破职业和地理区域之间的界限，成千上万的人会卖掉他们的房子，去拥抱一种在依靠互联网创造收入的同时周游世界的全新生活方式。这些人通过互联网赚取第一世界水平的收入，却选择生活在那些发展中国家物价水平的地方，他们被称作Digital Nomad（数字游民）。这种生活方式让他们彻底脱离了朝九晚五，办公室格挡和令人烦恼的通勤。”

如果你了解了“数字游民”这个概念，也有兴趣，有志于成为这样的数字游民，那么你现在所选择的《生物信息学》就是一个绝佳赛道。当然了，绝大部分小伙伴目前的技术实力谈不上独当一面，成为“数字游民”的路道阻且长，我这里有一个系统性提升自己能力的方案。

完成学徒作业，以markdown笔记的形式发到我邮箱，我会抽时间集中检查，挖掘其中足够优秀的小伙伴进行重点培养，给与更高级的学习资料或者个性化的学习指引，并且提供一定量的项目兼职测试一下你成为“数字游民”的潜力。

加油哦，我的邮箱是 jmzeng1314@163.com

文末友情推荐

做教学我们是认真的，如果你对我们的马拉松授课（直播一个月互动教学）有疑问，可以看完我们从2000多个提问互动交流里面精选的200个问答！ 2021第二期_生信入门班_微信群答疑整理，以及 2021第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信，你值得拥有下面的学习班：

github-actions[bot] commented 3 years ago

把单细胞表达量矩阵换一个单位 by 生信技能树

一般来说单细胞表达量矩阵都是以基因为单位，我们可以很容易走常规的降维聚类分群并且合理的进行生物学命名，比如我们对官方 pbmc3k 例子，跑标准代码：

library(Seurat)
# https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
## Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")

## Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", 
                           min.cells = 3, min.features = 200)
## Identification of mithocondrial genes
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

## Filtering cells following standard QC criteria.
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & 
                 percent.mt < 5)

## Normalizing the data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", 
                      scale.factor = 10000)

pbmc <- NormalizeData(pbmc)

## Identify the 2000 most highly variable genes
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

## In addition we scale the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), 
               verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)
pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")
table(pbmc$seurat_clusters)
# pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25,  logfc.threshold = 0.25, verbose = FALSE)
DimPlot(pbmc, reduction = "umap", group.by = 'seurat_clusters',
        label = TRUE, pt.size = 0.5) 
DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", 
                               "CD14", "FCER1A", "FCGR3A", 
                               "LYZ", "PPBP", "CD8A"),
        group.by = 'seurat_clusters')
## Assigning cell type identity to clusters
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T",
                     "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
 
pbmc$cluster_by_counts=Idents(pbmc)
table(pbmc$cluster_by_counts)

得到的结果如下所示：

降维聚类分群

根据我们的生物学背景，可以进行如下所示的命名：

生物学命名

这个时候的生物学背景因人而异哦，比如我这里对myeloid了解不多，所以就没有细分出来巨噬细胞，单核细胞，树突细胞等等。仅仅是举例而已。

首先可以替换成为转录因子活性

这里推荐一个 R包，DoRothEA on http://bioconductor.org/ and https://github.com/saezlab/dorothea.

获取包自带数据库

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("dorothea")
library(tidyverse)
## We read Dorothea Regulons for Human:
dorothea_regulon_human <- get(data("dorothea_hs", package = "dorothea"))

## We obtain the regulons based on interactions with confidence level A, B and C
regulon <- dorothea_regulon_human %>%
    dplyr::filter(confidence %in% c("A","B","C"))

使用数据库去分析pbmc

## We compute Viper Scores 
library(dorothea)
pbmc <- run_viper(pbmc, regulon,
                  options = list(method = "scale", minsize = 4, 
                                 eset.filter = FALSE, cores = 1, 
                                 verbose = FALSE))

重新进行降维聚类分群：

## We compute the Nearest Neighbours to perform cluster
DefaultAssay(object = pbmc) <- "dorothea"
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc, features = rownames(pbmc), verbose = FALSE)
pbmc <- FindNeighbors(pbmc, dims = 1:10, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)

pbmc <- RunUMAP(pbmc, dims = 1:10, umap.method = "uwot", metric = "cosine")
 
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()


library(gplots)
balloonplot(table(pbmc$dorothea_snn_res.0.5,
                  pbmc$cluster_by_counts))

同一个10x单细胞数据的两个矩阵各自的降维聚类分群是类似的，但是又不能完全一致：

这个时候的表达量矩阵并不是原始的每个基因在每个细胞里面的counts值了哦，而是每个转录因子在每个细胞的活性情况，我们也可以人工命名，

后面也可以对这样的矩阵（每个转录因子在每个细胞的活性情况），进行差异分析等等，参考；https://bioconductor.org/packages/release/data/experiment/vignettes/dorothea/inst/doc/single_cell_vignette.html

也可以可以替换成为gsea打分

单位就是各个基因集，比如MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据
C7: immunologic signatures: 免疫相关基因集合。

全新的矩阵是基因集在每个细胞的打分矩阵，然后可以再次聚类看看！

作为学徒作业

大家需要拿kegg的约300个通路，对这个pbmc3k单细胞表达量矩阵进行ssGSEA打分后得到矩阵，进行降维聚类分群！

写在文末

没有氛围，一个人吭哧吭哧的完成作业或多或少会遇到困难，如果有热火朝天的讨论，应该是会极大的促进咱们的积极性。所以有必要组建一个群聊，专门的针对于这个学徒作业100个练习题~

老规矩，群聊组建费用18.8元，一个简单的门槛隔绝那些不怀好意的广告营销号！ 前200名可以直接扫描（仍然需要18.8）群聊二维码进群，满200人后我们会统一收款！（每个人都是18.8 元，如果你不同意这个象征性收费，请不要进群哈！）

如果上面的二维码无法进群，说明已经满员了，需要我们生信技能树的官方拉群小助手帮忙拉群哦！！！（名额有限，先到先得！！！）这个时候请直接付款28元给小助手，就可以进群，或者你转发此推文到朋友圈然后截图给小助手，就可以仍然以18.8元进群！

一个简单的门槛，隔绝那些营销号！我们也会在群里共享学徒作业的资料，仅此而已，考虑清楚哦！

长按识别二维码

烦请备注姓名学校单位信息

会在群里共享收集整理好的全部文献

迈向数字游民之路

最近在朋友圈刷到了【一刻talks先见未来大会2020】的演讲视频，其中一个讲者很有意思、就是海尔集团海创汇合伙人，首席生态官——檀林，他的演讲：数字游民时代来临，你做好准备了吗？ 详细演讲视频介绍见：檀林：数字游民时代来临，你做好准备了吗？

让我接触到了“数字游民”这个词，其实最先提出“数字游民”构想的是一位名叫Tsugio Makimoto的日本高管。早在1997年，他便出版了一本名叫《Digital Nomad》的书：“未来的人类社会，高速的无线网络和强大的移动设备会打破职业和地理区域之间的界限，成千上万的人会卖掉他们的房子，去拥抱一种在依靠互联网创造收入的同时周游世界的全新生活方式。这些人通过互联网赚取第一世界水平的收入，却选择生活在那些发展中国家物价水平的地方，他们被称作Digital Nomad（数字游民）。这种生活方式让他们彻底脱离了朝九晚五，办公室格挡和令人烦恼的通勤。”

如果你了解了“数字游民”这个概念，也有兴趣，有志于成为这样的数字游民，那么你现在所选择的《生物信息学》就是一个绝佳赛道。当然了，绝大部分小伙伴目前的技术实力谈不上独当一面，成为“数字游民”的路道阻且长，我这里有一个系统性提升自己能力的方案。

完成学徒作业，以markdown笔记的形式发到我邮箱，我会抽时间集中检查，挖掘其中足够优秀的小伙伴进行重点培养，给与更高级的学习资料或者个性化的学习指引，并且提供一定量的项目兼职测试一下你成为“数字游民”的潜力。

加油哦，我的邮箱是 jmzeng1314@163.com

文末友情推荐

做教学我们是认真的，如果你对我们的马拉松授课（直播一个月互动教学）有疑问，可以看完我们从2000多个提问互动交流里面精选的200个问答！ 2021第二期_生信入门班_微信群答疑整理，以及 2021第二期_数据挖掘班_微信群答疑笔记

与十万人一起学生信，你值得拥有下面的学习班：

ixxmu commented 3 years ago

转录因子预测多了另一种软件预测，缺少 gsea

ixxmu / mp_duty

把单细胞表达量矩阵换一个单位 #1256

把单细胞表达量矩阵换一个单位 by 生信技能树

首先可以替换成为转录因子活性

获取包自带数据库

使用数据库去分析pbmc

也可以可以替换成为gsea打分

作为学徒作业

迈向数字游民之路

文末友情推荐

把单细胞表达量矩阵换一个单位 by 生信技能树

首先可以替换成为转录因子活性

获取包自带数据库

使用数据库去分析pbmc

也可以可以替换成为gsea打分

作为学徒作业

迈向数字游民之路

文末友情推荐