SCS【13】单细胞转录组之识别细胞对“基因集”的响应 (AUCell)_aucell 比较差异_桓峰基因的博客

软件包安装

if (!requireNamespace("BiocManager", quietly=TRUE))
  install.packages("BiocManager")
# To support paralell execution:
BiocManager::install(c("doMC", "doRNG","doSNOW"))
# For the main example:
BiocManager::install(c("mixtools", "SummarizedExperiment"))
# For the examples in the follow-up section of the tutorial:
BiocManager::install(c("DT", "plotly", "NMF", "d3heatmap", "shiny", "rbokeh",
                       "dynamicTreeCut","R2HTML","Rtsne", "zoo"))
AUCell所需要的输入数据分为两部分，矩阵数据和基因集数据。 
矩阵数据我们通过下载GEO上的数据集GSE60361下载，之后整理可以的矩阵： 
library(AUCell)
dir.create("AUCell_tutorial")
setwd("AUCell_tutorial")
if (!requireNamespace("GEOquery", quietly = TRUE)) BiocManager::install("GEOquery")
library(GEOquery)
# gse <- getGEO('GSE60361') # does not work, the matrix is in a suppl file
geoFile <- "GSE60361_C1-3005-Expression.txt.gz"
# download.file('https://ftp.ncbi.nlm.nih.gov/geo/series/GSE60nnn/GSE60361/suppl/GSE60361_C1-3005-Expression.txt.gz',
# destfile = geoFile)
library(data.table)
exprMatrix <- fread(geoFile, sep = "\t")
geneNames <- unname(unlist(exprMatrix[, 1, with = FALSE]))
exprMatrix <- as.matrix(exprMatrix[, -1, with = FALSE])
rownames(exprMatrix) <- geneNames
exprMatrix <- exprMatrix[unique(rownames(exprMatrix)), ]
dim(exprMatrix)
exprMatrix[1:5, 1:4]
# Remove file(s) downloaded:
file.remove(geoFile)
# Convert to sparse
library(Matrix)
exprMatrix <- as(exprMatrix, "dgCMatrix")
# Save for future use
mouseBrainExprMatrix <- exprMatrix
save(mouseBrainExprMatrix, file = "exprMatrix_MouseBrain.RData")
基因集是AUCell软件包自带的直接调取利用getGmt读取即可： 
library(AUCell)
library(GSEABase)
gmtFile <- paste(file.path(system.file("examples", package = "AUCell")), "geneSignatures.gmt",
    sep = "/")
geneSets <- getGmt(gmtFile)
geneSets <- subsetGeneSets(geneSets, rownames(exprMatrix))
cbind(nGenes(geneSets))
geneSets <- setGeneSetNames(geneSets, newNames = paste(names(geneSets), " (", nGenes(geneSets),
    "g)", sep = ""))
# Random
set.seed(321)
extraGeneSets <- c(GeneSet(sample(rownames(exprMatrix), 50), setName = "Random (50g)"),
    GeneSet(sample(rownames(exprMatrix), 500), setName = "Random (500g)"))
countsPerGene <- apply(exprMatrix, 1, function(x) sum(x > 0))
# Housekeeping-like
extraGeneSets <- c(extraGeneSets, GeneSet(sample(names(countsPerGene)[which(countsPerGene >
    quantile(countsPerGene, probs = 0.95))], 100), setName = "HK-like (100g)"))
geneSets <- GeneSetCollection(c(geneSets, extraGeneSets))
names(geneSets)
计算gene signatures分值
 
为每个细胞建立基因表达排序，然后进行计算基因标记的富集(AUC)： 
cells_AUC <- AUCell_run(exprMatrix, geneSets)
save(cells_AUC, file = "cells_AUC.RData")
cells_rankings <- AUCell_buildRankings(exprMatrix, plotStats = FALSE, splitByBlocks = TRUE)
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings)
cells_rankings
确定具有给定基因特征或活性基因集的细胞
 
AUC表示在signature中表达的基因的比例，以及与细胞内其他基因的相对表达值。我们可以使用此属性根据基因集的表达来探索数据集中存在的细胞种群。 
set.seed(123)
par(mfrow = c(3, 3))
cells_assignment <- AUCell_exploreThresholds(cells_AUC, plotHist = TRUE, assign = TRUE)
绘制特定基因集的AUC直方图，并设置新的阈值: 
cells_assignment$Oligodendrocyte_Cahoy$aucThr$thresholds
cells_assignment$Oligodendrocyte_Cahoy$aucThr$selected
oligodencrocytesAssigned <- cells_assignment$Oligodendrocyte_Cahoy$assignment
length(oligodencrocytesAssigned)
head(oligodencrocytesAssigned)
geneSetName <- rownames(cells_AUC)[grep("Oligodendrocyte_Cahoy", rownames(cells_AUC))]
AUCell_plotHist(cells_AUC[geneSetName, ], aucThr = 0.25)
abline(v = 0.25)
为这个新的阈值分配细胞: 
newSelectedCells <- names(which(getAUC(cells_AUC)[geneSetName, ] > 0.08))
length(newSelectedCells)
head(newSelectedCells)
探索细胞分配(表格和热图)，细胞赋值存储在$assignment中。提取所有基因集的细胞，并将其转化为一个表格，并绘制为直方图: 
cellsAssigned <- lapply(cells_assignment, function(x) x$assignment)
assignmentTable <- reshape2::melt(cellsAssigned, value.name = "cell")
colnames(assignmentTable)[2] <- "geneSet"
head(assignmentTable)
assignmentMat <- table(assignmentTable[, "geneSet"], assignmentTable[, "cell"])
assignmentMat[, 1:2]
set.seed(123)
miniAssigMat <- assignmentMat[, sample(1:ncol(assignmentMat), 100)]
library(NMF)
aheatmap(miniAssigMat, scale = "none", color = "black", legend = FALSE)
library(DT)
datatable(assignmentTable, options = list(pageLength = 10), filter = "top")
根据signatrue评分探索细胞/簇
 
AUC评分还可以用于探索以前聚类分析的输出(反之亦然)。 
使用获得的不同signature的AUC为t-SNE上色，该t-SNE基于之前运行的整个表达矩阵(即不是5000个随机基因)。 
load(paste(file.path(system.file("examples", package = "AUCell")), "cellsTsne.RData",
    sep = "/"))
cellsTsne <- cellsTsne$Y
plot(cellsTsne, pch = 16, cex = 0.3)
这个t-SNE是用下面的代码创建的(需要一段时间来运行): 
sumByGene <- apply(mouseBrainExprMatrix, 1, sum)
exprMatSubset <- mouseBrainExprMatrix[which(sumByGene > 0), ]
logMatrix <- log2(exprMatSubset + 1)
logMatrix <- as.matrix(logMatrix)
library(Rtsne)
set.seed(123)
cellsTsne <- Rtsne(t(logMatrix))  ##运行时间长
rownames(cellsTsne$Y) <- colnames(logMatrix)
colnames(cellsTsne$Y) <- c("tsne1", "tsne2")
save(cellsTsne, file = "cellsTsne.RData")
load("cellsTsne.RData")
t-SNE可以根据AUC评分进行着色。为了突出显示根据特征更可能属于该细胞类型的细胞簇，我们将把细胞分为通过分配阈值的细胞(用粉色-红色表示)和没有通过分配阈值的细胞(用黑色-蓝色表示)。 
当然，在解释时也应该记住这些签名的来源(例如，这些签名是从大体积RNA-seq分析中获得的)。此外，只在5000个基因上运行本教程可能会带来一些额外的噪音。 
对于本例，我们使用了自动分配的阈值: 
selectedThresholds <- getThresholdSelected(cells_assignment)
par(mfrow = c(3, 3))  # Splits the plot into two rows and three columns
for (geneSetName in names(selectedThresholds)) {
    nBreaks <- 5# Number of levels in the color palettes
    # Color palette for the cells that do not pass the threshold
    colorPal_Neg <- (grDevices::colorRampPalette(c("black", "blue", "skyblue")))(nBreaks)
    # Color palette for the cells that pass the threshold
    colorPal_Pos <- (grDevices::colorRampPalette(c("pink", "magenta", "red")))(nBreaks)
    # Split cells according to their AUC value for the gene set
    passThreshold <- getAUC(cells_AUC)[geneSetName, ] > selectedThresholds[geneSetName]
    if (sum(passThreshold) > 0) {
        aucSplit <- split(getAUC(cells_AUC)[geneSetName, ], passThreshold)
        # Assign cell color
        cellColor <- c(setNames(colorPal_Neg[cut(aucSplit[[1]], breaks = nBreaks)],
            names(aucSplit[[1]])), setNames(colorPal_Pos[cut(aucSplit[[2]], breaks = nBreaks)],
            names(aucSplit[[2]])))
        # Plot
        plot(cellsTsne$Y, main = geneSetName, sub = "Pink/red cells pass the threshold",
            col = cellColor[rownames(cellsTsne$Y)], pch = 16)
这种图也可以用来探索分配阈值。例如: 
selectedThresholds[2] <- 0.25
par(mfrow = c(2, 3))
AUCell_plotTSNE(tSNE = cellsTsne$Y, exprMat = exprMatrix, cellsAUC = cells_AUC[1:2,
    ], thresholds = selectedThresholds)
与均值比较
 
评估基因特征表达最直观的方法之一是计算平均表达。然而，计算原始计数上的基因签名的平均值将导致读取次数更多的细胞(无论是技术上还是生物学上的原因)拥有更高的大多数签名的平均值。为了减少这种影响，需要某种类型的库大小规范化。由于AUCell单独评估每个单元上的签名检测，它自动补偿库的大小，并使其应用程序更容易，独立于数据单位和应用的其他规范化(例如原始计数、TPM、UMI计数)。 
## Comparison with mean
logMat <- log2(exprMatrix + 2)
meanByGs <- t(sapply(geneSets, function(gs) colMeans(logMat[geneIds(gs), ])))
rownames(meanByGs) <- names(geneSets)
colorPal <- (grDevices::colorRampPalette(c("black", "red")))(nBreaks)
par(mfrow = c(3, 3))
for (geneSetName in names(geneSets)) {
    cellColor <- setNames(colorPal[cut(meanByGs[geneSetName, ], breaks = nBreaks)],
        names(meanByGs[geneSetName, ]))
    plot(cellsTsne$Y, main = geneSetName, axes = FALSE, xlab = "", ylab = "", sub = "Expression mean",
        col = cellColor[rownames(cellsTsne$Y)], pch = 16)
绘制tSNE降维的AUC： 
AUCell_plotTSNE(tSNE = cellsTsne$Y, exprMat = exprMatrix, plots = "AUC", cellsAUC = cells_AUC[1:9,
    ], thresholds = selectedThresholds)
检测到的基因数量 
nGenesPerCell <- apply(exprMatrix, 2, function(x) sum(x > 0))
colorPal <- grDevices::colorRampPalette(c("darkgreen", "yellow", "red"))
cellColorNgenes <- setNames(adjustcolor(colorPal(10), alpha.f = 0.8)[as.numeric(cut(nGenesPerCell,
    breaks = 10, right = FALSE, include.lowest = TRUE))], names(nGenesPerCell))
plot(cellsTsne$Y, axes = FALSE, xlab = "", ylab = "", sub = "Number of detected genes",
    col = cellColorNgenes[rownames(cellsTsne$Y)], pch = 16)
How reliable is AUCell? (Confusion matrix)
 
当使用AUCell识别单元格类型时，所识别单元格的可靠性将取决于使用的签名，以及数据集中单元格类型之间的差异。在任何情况下，AUCell对噪声都相当稳健，无论是在数据上还是在基因集上。 
这是用AUCell(本例)识别的细胞与最初发表的细胞类型的比较: 
# 'Real' cell type (e.g. provided in the publication)
cellLabels <- paste(file.path(system.file("examples", package = "AUCell")), "mouseBrain_cellLabels.tsv",
    sep = "/")
cellLabels <- read.table(cellLabels, row.names = 1, header = TRUE, sep = "\t")
# Confusion matrix:
cellTypeNames <- unique(cellLabels[, "level1class"])
confMatrix <- t(sapply(cells_assignment[c(1, 4, 2, 5, 3, 6:9)], function(x) table(cellLabels[x$assignment,
    ])[cellTypeNames]))
colnames(confMatrix) <- cellTypeNames
confMatrix[which(is.na(confMatrix), arr.ind = TRUE)] <- 0
confMatrix
我们这期主要介绍单细胞转录组数据中识别细胞对“基因集”的响应 (AUCell)。目前单细胞测序的费用也在降低，单细胞系列可算是目前的测序神器，有这方面需求的老师，联系桓峰基因，提供最高端的科研服务！ 
桓峰基因，铸造成功的您！ 
未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程， 
敬请期待！！ 
有想进生信交流群的老师可以扫最后一个二维码加微信，备注“单位+姓名+目的”，有些想发广告的就免打扰吧，还得费力气把你踢出去！ 
References:
 
 Aibar et al. (2017) SCENIC: single-cell regulatory network inference and clustering. Nature Methods. doi: 10.1038/nmeth.4463
 
 Cahoy, J.D., et al. (2008). A Transcriptome Database for Astrocytes, Neurons, and Oligodendrocytes: A New Resource for Understanding Brain Development and Function. J. Neurosci. 28, 264–278. doi: 10.1523/JNEUROSCI.4178-07.2008
 
 Lein, E.S., et al. (2007). Genome-wide atlas of gene expression in the adult mouse brain. Nature 445, 168–176. doi: 10.1038/nature05453
 
 Lavin, Y., et al. (2014) Tissue-Resident Macrophage Enhancer Landscapes Are Shaped by the Local Microenvironment. Cell 159, 1312–1326. doi: 10.1016/j.cell.2014.11.018
  
				更多精彩请至我的公众号《KS科研分享与服务》
第一次接触基因评分是在一篇文章中，也不知道这样的叫法对不对，作者选定了几个炎症基因，利用seurat包的一个打分函数AddModuleScore，依据基因的平均表达水平进行分析，最后得到的score称为炎症分数，其实这样的叫法有点欠缺，但是有这样的做法。（演示数据没有意义）
选择基因进行计算。
DefaultAssay(immune) <- "RNA"
cd_features <- list(c(
  'TNF',
  'CCL2',
				第一次接触基因评分是在一篇文章中，也不知道这样的叫法对不对，作者选定了几个炎症基因，利用seurat包的一个打分函数AddModuleScore，依据基因的平均表达水平进行分析，最后得到的score称为炎症分数，其实这样的叫法有点欠缺，但是有这样的做法。（演示数据没有意义）
选择基因进行计算。
DefaultAssay(immune) <- “RNA”
cd_features <- list(c(
‘TNF’,
‘CCL2’,
‘CCL3’,
‘CCL4’,
‘CXCL10’,
‘S100A8’
				基因签名发现
 基于乳腺癌亚型患者的microRNA，mRNA定量数据和其他临床数据，我们希望通过机器学习建立模型以预测5年生存率和肿瘤分期。 如果不能正确预测这些标签，我们希望通过聚类进行数据挖掘，以找出是否可以将数据分类为其他子类以及区分这些子类的基因。 此外，我们想找出那些基因参与的生物学途径，并对微RNA-mRNA对的特征做出假设。
 - Writen report can be found here: (https://github.com/timothyfisherphd/Gene_Signature_Discovery/blob/main/Gene_Signature_Discovery.pdf)
打开终端或命令窗口以执行以下代码：./run.sh
 使用基因列表的输出，然后通过/ Pathway_Analysis运行基因名称。 打开Python3并在./GO_
				相信做过肿瘤单细胞的小伙伴对这个分析并不陌生，如果多读几篇文献，就能在CNS以及大子刊上面看到这个分析。
非负矩阵分解(Nonnegative Matrix Factorization,NMF)是在矩阵中所有元素均为非负数约束条件之下的矩阵分解方法。
基本思想：给定一个非负矩阵V, NMF能够找到一个非负矩阵W和一个非负矩阵H, 使得矩阵W和H的乘积近似等于矩阵V中的值。
放在我们单细胞转录组的场景下，就是需要将一个基因×细胞的表达矩阵(V)，分解成基因×表达程序(W)，与表达程序×细胞(H)两个矩阵的
				R读取gmt文件的五种方式
kegg_geneset1 <- qusage::read.gmt("c2.cp.kegg.v7.0.symbols.gmt") #返回的是list
kegg_geneset2 <- clusterProfiler::read.gmt("c2.cp.kegg.v7.0.symbols.gmt") #返回的是表格
kegg_geneset3 <- GSA::GSA.read.gmt("c2.cp.kegg.v7.0.symbols.gmt") #"GSA.gene
				SCS（Subcarrier Spacing，子载波间隔）是5G信号的基本参数之一，用于表示子载波间隔大小，通常用kHz作为单位。SCS不同取值对应着不同的时域资源格子大小，即用于时域调度的最小时间单元（TTI，Transmission Time Interval）。在5G NR（New Radio）中，SCS的取值有15 kHz、30 kHz、60 kHz、120 kHz、240 kHz五种，它们与时域的对应关系如下：
- 当SCS取值为15 kHz时，时域资源格子大小为0.5 ms。
- 当SCS取值为30 kHz时，时域资源格子大小为1 ms。
- 当SCS取值为60 kHz时，时域资源格子大小为2 ms。
- 当SCS取值为120 kHz时，时域资源格子大小为4 ms。
- 当SCS取值为240 kHz时，时域资源格子大小为8 ms。
从上述对应关系可以看出，SCS取值越大，时域资源格子越大，也就意味着传输速率越快，但灵活性更差。因此，在5G系统中需要根据具体的应用场景选择不同的SCS取值。例如，对于高速移动的用户，建议使用SCS取值较大的情况下，因为较大的SCS可以提高传输速率，适应高速移动的需求。而在室内、短距离通信等场景下，可以采用SCS取值较小的情况，以支持更高的信道容量和更好的低延迟性能。