vst 算法R语言手工实现 | Seurat4 筛选高变基因的算法

1. vst算法描述

(1)为什么需要矫正

在这里插入图片描述
image source: https://ouyanglab.com/singlecell/basic.html

In this panel, we observe that there is a very strong positive relationship between a gene’s average expression and its observed variance. In other words, highly expressed genes have high variances associated with them and vice versa. This phenomenon is often referred to as heteroskedasticity within the data and must be corrected before proceeding with any downstream analysis steps.
高表达的基因其变异也高。这被叫做数据内的 异方差性,必须矫正后才能进一步下游分析。

In short, the dependance of a gene’s expression with its observed variance is what needs to be corrected prior to HVG selection [1].
简而言之,基因表达与其观察到的变异的依赖性是在HVG选择之前需要纠正的。

The Solution:

A very common way to correct this problem is by applying a Variance Stabilizing Transformation (VST) to the data, and we see in the right panel of the above figure that once VST is applied, the relationship of observed variance at any given level of average expression for a gene has been removed/standardized.

(2)Seurat4 R文档

> ?FindVariableFeatures
vst:

  • First, fits a line to the relationship of log(variance) and log(mean) using local polynomial regression (loess).
  • Then standardizes the feature values using the observed mean and expected variance (given by the fitted line).
  • Feature variance is then calculated on the standardized values after clipping to a maximum (see clip.max parameter).

vst steps: 目的是在var~mean曲线中,不同mean值区域都能挑选var较大的基因

  1. 使用局部多项式拟合(loess) log(variance) 和log(mean) 平滑曲线模型
  2. 获取模型计算的值作为y=var.exp值
  3. 截取最大之后,var.standarlized = get variance after feature standardization:
    (每个基因 - mean)/sd 后 取var(). 注意sd=sqrt(var.exp)
  4. 按照 var.standarlized 降序排列,取前n(比如2000)个基因作为高变基因。

2. 加载数据及Seurat vst 结果

使用pbmc 3k数据,走标准Seurat4,选取top 2000 HVG [2]。

library(Seurat)
library(ggplot2)
library(dplyr)pbmc=readRDS("d:\\code_R\\filtered_gene_bc_matrices\\pbmc3k.final.Rds")
DimPlot(pbmc, label=T)# pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)# 0. 获取Seurat包计算的HVG ----
top10 <- head(VariableFeatures(pbmc), 10)
p1=VariableFeaturePlot(pbmc); p1
LabelPoints(plot = p1, points = top10, repel = TRUE)genelist1=VariableFeatures(pbmc)
length(genelist1)
head(genelist1)head( pbmc@assays$RNA@meta.features )

3. 手工 HVG

(1) LOESS fit y~x

log(variance) ~ log(mean)

# use raw data: vst 的直接输入的是 counts,所以直接算的 行平均数,作为基因表达值
counts = pbmc@assays$RNA@counts# 计算每行均值
hvf.info <- data.frame(mean = rowMeans(x = counts, na.rm = T))# 计算每行方差
hvf.info$variance = apply(counts, 1, function(x){var(x, na.rm = T)
})head(hvf.info)
dim(hvf.info)if(0){#pdf( paste0(outputRoot, keyword, "_01_4B.HVG.pdf"), width=5.5, height =4.8)plot(hvf.info$mean, hvf.info$variance, pch=19, cex=0.3)plot(log10(hvf.info$mean), hvf.info$variance, pch=19, cex=0.3)plot(log10(hvf.info$mean), log10(hvf.info$variance), pch=19, cex=0.3) #looks good#dev.off()
}# 通过loess拟合,计算期望方差和标准化的方差
hvf.info$variance.expected <- 0
hvf.info$variance.standardized <- 0not.const <- (hvf.info$variance >0) & (!is.na(hvf.info$variance)) & (!is.na(hvf.info$mean))
table(not.const)
#TRUE 
#13714# 拟合 y~x
loess.span=0.3 #default in Seurat4
fit <- loess(formula = log10(x = variance) ~ log10(x = mean),data = hvf.info[not.const, ],span = loess.span
)dim(hvf.info[not.const, ]) #13714     4
#str(fit)

(2). 获取模型给出的期望y值

# 期望y值:使用模型计算的值
hvf.info$variance.expected[not.const] <- 10 ^ fit$fitted

(3). 截取极端值后,计算标准化后的方差

#clip.max == 'auto',则自动设置为 列数(细胞数)的开方
clip.max <- sqrt(x = ncol(x = counts))
clip.max #51.9# 计算feature标准化( (counts - mean)/sd )后的方差,注意sd=sqrt(var)
# 求方差前截取极大值
hvf.info$variance.standardized[not.const]=(function(){result=c()for(i in 1:nrow(counts)){row.var=NAif(not.const[i]){# clip to a maximumrow.counts.std = (counts[i, ] - hvf.info$mean[i]) / sqrt(hvf.info$variance.expected[i])row.counts.std[row.counts.std>clip.max]=clip.max# 计算方差row.var=var( row.counts.std, na.rm = T )}# 返回结果result=c(result, row.var)}return(result)
})() #2min

(4) 第(3)步主要参考

来源是Seurat4 函数FindVariableFeatures.default 中:

hvf.info$variance.standardized <- SparseRowVarStd(mat = object,mu = hvf.info$mean,sd = sqrt(hvf.info$variance.expected),vmax = clip.max,display_progress = verbose)

查找SparseRowVarStd的定义:

$ find .| xargs grep -in "SparseRowVarStd" --color=auto
grep: .: Is a directory
./data_manipulation.cpp:305:NumericVector SparseRowVarStd(Eigen::SparseMatrix<double> mat,
./data_manipulation.h:40:NumericVector SparseRowVarStd(Eigen::SparseMatrix<double> mat,
./RcppExports.cpp:185:// SparseRowVarStd
./RcppExports.cpp:186:NumericVector SparseRowVarStd(Eigen::SparseMatrix<double> mat, NumericVector mu, NumericVector sd, double vmax, bool display_progress);
./RcppExports.cpp:187:RcppExport SEXP _Seurat_SparseRowVarStd(SEXP matSEXP, SEXP muSEXP, SEXP sdSEXP, SEXP vmaxSEXP, SEXP display_progressSEXP) {
./RcppExports.cpp:195:    rcpp_result_gen = Rcpp::wrap(SparseRowVarStd(mat, mu, sd, vmax, display_progress));
./RcppExports.cpp:421:    {"_Seurat_SparseRowVarStd", (DL_FUNC) &_Seurat_SparseRowVarStd, 5},

发现是Seurat4 的 c++函数,定义在 seurat-4.1.0/src/data_manipulation.cpp:301

/* standardize matrix rows using given mean and standard deviation,clip values larger than vmax to vmax,then return variance for each row */
// [[Rcpp::export(rng = false)]]
NumericVector SparseRowVarStd(Eigen::SparseMatrix<double> mat,NumericVector mu,NumericVector sd,double vmax,bool display_progress){if(display_progress == true){Rcpp::Rcerr << "Calculating feature variances of standardized and clipped values" << std::endl;}mat = mat.transpose();NumericVector allVars(mat.cols());Progress p(mat.outerSize(), display_progress);for (int k=0; k<mat.outerSize(); ++k){p.increment();if (sd[k] == 0) continue;double colSum = 0;int nZero = mat.rows();for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it){nZero -= 1;colSum += pow(std::min(vmax, (it.value() - mu[k]) / sd[k]), 2);}colSum += pow((0 - mu[k]) / sd[k], 2) * nZero;allVars[k] = colSum / (mat.rows() - 1);}return(allVars);
}

结合R的上下文,大致能猜出来c++代码啥意思。
想看懂细节,则需要更多c++知识储备:

  • cpp最顶部的头文件:
#include <RcppEigen.h>
#include <progress.hpp>
#include <cmath>
#include <unordered_map>
#include <fstream>
#include <string>
#include <Rinternals.h>using namespace Rcpp;
  • Eigen::SparseMatrix<double> mat:泛型;矩阵::系数矩阵
  • 稀疏矩阵的方法:mat.transpose(), mat.cols(), mat.outerSize(),
  • for (Eigen::SparseMatrix<double>::InnerIterator it(mat,k); it; ++it): 集合的迭代器
  • Rcpp 数据类型:NumericVector,

4. 结果比较

(1) 结果检查1:基因列表一致

手工计算的和Seurat4的HVG gene list结果完全一致。

# Check 1: HVG gene list
hvf.info=hvf.info[order(-hvf.info$variance.standardized),]
head(hvf.info)
tail(hvf.info)
#top.features=head( rownames(hvf.info), n=250)
top.features_2=head( rownames(hvf.info), n=2000)
length(top.features_2)
head(top.features_2)
setdiff(top.features_2, genelist1)
setdiff(genelist1, top.features_2)
#
if(0){# Check: gene and their parampbmc@assays$RNA@meta.features[c(setdiff(genelist1, top.features_2)),]hvf.info[c(setdiff(genelist1, top.features_2)),]#pbmc@assays$RNA@meta.features[c(setdiff(top.features_2, genelist1)),]hvf.info[c(setdiff(top.features_2, genelist1)),]
}

(2) 结果检查2:基因参数一致

手工计算的和Seurat4的HVG gene 参数完全一致。

# check2: HVG and its params
# 1.
dim(pbmc@assays$RNA@meta.features)
head( pbmc@assays$RNA@meta.features )
#                 vst.mean vst.variance vst.variance.expected vst.variance.standardized vst.variable
#AL627309.1    0.003333333  0.003323453           0.003575582                 0.9294859        FALSE
#AP006222.2    0.001111111  0.001110288           0.001112798                 0.9977442        FALSE
#RP11-206L10.2 0.001851852  0.001849107           0.001921811                 0.9621691        FALSE
#RP11-206L10.9 0.001111111  0.001110288           0.001112798                 0.9977442        FALSE
#LINC00115     0.006666667  0.006624676           0.007342308                 0.9022607        FALSE
#NOC2L         0.106666667  0.158310485           0.203482316                 0.7780061        FALSE# 2.
dim(hvf.info)
hvf.info=hvf.info[rownames(pbmc@assays$RNA@meta.features),]
head(hvf.info)
#                     mean    variance variance.expected variance.standardized
#AL627309.1    0.003333333 0.003323453       0.003575582             0.9294859
#AP006222.2    0.001111111 0.001110288       0.001112798             0.9977442
#RP11-206L10.2 0.001851852 0.001849107       0.001921811             0.9621691
#RP11-206L10.9 0.001111111 0.001110288       0.001112798             0.9977442
#LINC00115     0.006666667 0.006624676       0.007342308             0.9022607
#NOC2L         0.106666667 0.158310485       0.203482316             0.7780061

比较高变基因的参数:

> hvf.info[genelist1|> head(), ]mean    variance variance.expected variance.standardized
PPBP    0.2451852    9.577506         0.5888573             11.171153
S100A9  6.0466667  278.681037        34.8969051              7.985838
IGLL5   0.2792593    8.894938         0.6929479              7.964379
LYZ    10.2466667  564.108825        70.8492711              7.962098
GNLY    1.5740741   45.239046         6.0378423              7.492585
FTL    27.6674074 2008.688897       278.9968379              7.199683
> pbmc@assays$RNA@meta.features[genelist1|> head(), ]vst.mean vst.variance vst.variance.expected vst.variance.standardized vst.variable
PPBP    0.2451852     9.577506             0.5888573                 11.172765         TRUE
S100A9  6.0466667   278.681037            34.8969051                  7.985838         TRUE
IGLL5   0.2792593     8.894938             0.6929479                  7.965360         TRUE
LYZ    10.2466667   564.108825            70.8492711                  7.962098         TRUE
GNLY    1.5740741    45.239046             6.0378423                  7.492585         TRUE
FTL    27.6674074  2008.688897           278.9968379                  7.199683         TRUE> table(abs(hvf.info[genelist1, 4] - pbmc@assays$RNA@meta.features[genelist1, 4])<0.005)
FALSE  TRUE 3  1997# 差别不大,差异的绝对值大于0.005的共三个:
> keep2=abs(hvf.info[genelist1, 4] - pbmc@assays$RNA@meta.features[genelist1, 4])>0.005
> table(keep2)
keep2
FALSE  TRUE 1997     3 > hvf.info[genelist1, ][keep2, ]mean  variance variance.expected variance.standardized
IGJ      0.16777778 16.822896        0.36540081              3.481455
SLC48A1  0.03370370  0.871409        0.04618194              2.215032
NAPA-AS1 0.02962963  1.050622        0.03932270              1.274513
> pbmc@assays$RNA@meta.features[genelist1, ][keep2, ]vst.mean vst.variance vst.variance.expected vst.variance.standardized vst.variable
IGJ      0.16777778    16.822896            0.36540081                  3.498952         TRUE
SLC48A1  0.03370370     0.871409            0.04618194                  2.220942         TRUE
NAPA-AS1 0.02962963     1.050622            0.03932270                  1.280866         TRUE

最后一列略有区别:

#
all( abs(hvf.info[,1] - pbmc@assays$RNA@meta.features[,1]) < 1e-6) #T
all( abs(hvf.info[,2] - pbmc@assays$RNA@meta.features[,2]) < 1e-6) #T
all( abs(hvf.info[,3] - pbmc@assays$RNA@meta.features[,3]) < 1e-6) #Ttable( abs(hvf.info[,4] - pbmc@assays$RNA@meta.features[,4]) < 1e-6) #not all T
table( abs(hvf.info[,4] - pbmc@assays$RNA@meta.features[,4]) < 0.01)
#FALSE  TRUE 
#    1 13713keep = abs(hvf.info[,4] - pbmc@assays$RNA@meta.features[,4]) > 1e-2> table(keep)
keep
FALSE  TRUE 
13713     1 
> hvf.info[keep, ]mean variance variance.expected variance.standardized
IGJ 0.1677778  16.8229         0.3654008              3.481455
> pbmc@assays$RNA@meta.features[keep, ]vst.mean vst.variance vst.variance.expected vst.variance.standardized vst.variable
IGJ 0.1677778      16.8229             0.3654008                  3.498952         TRUE

(3) 绘图比较

(a) (a) var~avg with top 2000 HVG selected by vst;
(b) std.var ~ avg with top 2000 HVG selected by vst;
( c) same as (b), but draw by Seurat functions.

# plot1
plot(log10(hvf.info$mean), log10(hvf.info$variance), pch=19, cex=0.3, main="vst manully #1", mgp=c(2,1,0))
points(log10(hvf.info[top.features_2, ]$mean), log10(hvf.info[top.features_2, ]$variance), pch=19, cex=0.3, col="red")# plot2
plot(log10(hvf.info$mean), hvf.info$variance.standardized, pch=19, cex=0.3, main="vst manully #2", mgp=c(2,1,0))
points(log10(hvf.info[top.features_2,]$mean), (hvf.info[top.features_2,]$variance.standardized), pch=19, cex=0.3, col="red")# plot3: Seurat
LabelPoints(plot = p1, points = top10, repel = TRUE)
#

Ref:

  • [1] https://medium.com/byte-sized-machine-learning/selection-of-highly-variable-genes-hvgs-in-scrna-seq-647c8eee3845
  • [2] https://zhuanlan.zhihu.com/p/479549742

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/bicheng/48166.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

OpenAI 推出 GPT-4o mini,一种更小、更便宜的人工智能模型

OpenAI 最近推出了新型人工智能模型 GPT-4o mini&#xff0c;以其较小体积和低成本受到关注。这款模型在文本和视觉推理任务上性能优越&#xff0c;且比现有小型模型更快、更经济。GPT-4o mini 已向开发者和消费者发布&#xff0c;企业用户将在下周获得访问权限。 喜好儿网 在…

入坑树莓派(2)——树莓派4B与手机蓝牙通信

入坑树莓派(2)——树莓派4B与手机蓝牙通信 1、引言 在入坑树莓派(1)中已经搞掂了可视化问题。现在继续开展下一步,尝试与手机通信,一开始是想弄wifi连接的,但发现基于wifi的APP比较难弄,为了降低开发的难度,又因为树莓派板子自带蓝牙模块,所以直接选用蓝牙连接手机…

LabVIEW多串口通信

随着现代工业控制对数据采集和处理效率的要求不断提升&#xff0c;传统的单串口通信已无法满足多通道数据传输与大规模数据存取的需求。开发一种基于LabVIEW的多串口通信及数据存储系统&#xff0c;以提升数据处理速度和存储效率&#xff0c;保障生产线的稳定运行显得尤为重要。…

数据结构之树的存储结构详解与示例(C/C++)

文章目录 树的存储结构1. 顺序存储结构2. 链式存储结构结论 树&#xff08;Tree&#xff09;是一种非常常见的数据结构&#xff0c;它模拟了一种层级或分支结构。树由节点&#xff08;或称为顶点&#xff09;组成&#xff0c;每个节点包含一个值&#xff0c;并且可能有多个子节…

SpringDoc2问题汇总

在项目中尝试使用SpringDoc进行文档生成&#xff0c;在使用过程中遇到一系列的问题加以记录. 1.引入依赖 只是单纯的使用SpringDoc的话不需要引入一些乱七八糟的依赖&#xff0c;如今各种增强和拓展依赖层出不穷&#xff0c;但是随着这些依赖的出现带来的不仅是增强&#xff0…

在学习使用LabVIEW的过程中,需要注意哪些问题?

在学习使用LabVIEW的过程中&#xff0c;需要注意以下问题&#xff1a; 1. 基础知识 图形化编程思维&#xff1a; LabVIEW采用图形化编程方式&#xff0c;与传统的文本编程语言有很大不同&#xff0c;需要适应这种新的编程思维方式。数据流概念&#xff1a; 理解LabVIEW的核心数…

服务器借助笔记本热点WIFI上网

一、同一局域网环境 1、当前环境&#xff0c;已有交换机组网环境&#xff0c;服务器已配置IP信息。 设备ip服务器125.10.100.12交换机125.10.100.0/24笔记本125.10.100.39 2、拓扑图 #mermaid-svg-D4moqMym9i0eeRBm {font-family:"trebuchet ms",verdana,arial,sa…

AFAC2024-基于保险条款的问答 比赛日记 llamafactory qwen npu 910B1

AFAC2024: 基于保险条款的问答挑战——我的实战日记 概述 在最近的AFAC2024竞赛中&#xff0c;我参与了基于保险条款的问答赛道。这是一次深度学习与自然语言处理的实战演练&#xff0c;旨在提升模型在复杂保险文本理解与问答生成方面的能力。本文将分享我的参赛过程&#xf…

【权威发布】2024年生物技术与医学国际会议(IACBM 2024)

2024年生物技术与医学国际会议 2024 International Conference on Biotechnology and Medicine 【1】会议简介 2024年生物技术与医学国际会议旨在为全球生物技术与医学领域的专家学者提供一个交流最新研究成果、分享技术进展和探讨未来发展方向的平台。会议旨在加强国际间的学术…

阿里云 https证书部署

一.申请证书 二.查看状态 查看状态&#xff0c;已签发是完成了申请证书 三.部署 我在nginx服务器上部署 具体操作链接:阿里云文档 修改前 修改后 四.重启ngnix 五.验证是否成功 在浏览器输入域名查看

【JavaScript 算法】最长公共子序列:字符串问题的经典解法

&#x1f525; 个人主页&#xff1a;空白诗 文章目录 一、算法原理状态转移方程初始条件 二、算法实现注释说明&#xff1a; 三、应用场景四、总结 最长公共子序列&#xff08;Longest Common Subsequence&#xff0c;LCS&#xff09;是字符串处理中的经典问题。给定两个字符串…

ETL电商项目总结

ETL电商项目总结 ETL电商业务简介及各数据表关系 业务背景 ​ 本案例围绕某个互联网小型电商的订单业务来开发。某电商公司&#xff0c;每天都有一些的用户会在线上采购商品&#xff0c;该电商公司想通过数据分析&#xff0c;查看每一天的电商经营情况。例如&#xff1a;电商…

通信协议_C#实现CAN通信

CAN协议 CAN&#xff08;Controller Area Network&#xff09;即控制器局域网络。特点&#xff1a; 多主网络&#xff1a;网络上的任何节点都可以主动发送数据&#xff0c;不需要一个固定的主节点。双绞线&#xff1a;使用双绞线作为通信介质&#xff0c;支持较远的通信距离。…

时序数据库如何选型?详细指标总结!

工业物联网场景&#xff0c;如何判断什么才是好的时序数据库&#xff1f; 工业物联网将机器设备、控制系统与信息系统、业务过程连接起来&#xff0c;利用海量数据进行分析决策&#xff0c;是智能制造的基础设施&#xff0c;并影响整个工业价值链。工业物联网机器设备感知形成了…

C++那些事之依赖注入

C那些事之依赖注入 最近星球里面有个小伙伴让更新一下依赖注入&#xff0c;于是写出了这篇文章&#xff0c;来从实际的例子讲解&#xff0c;本文会讲解一些原理与实现&#xff0c;完整的实现代码懒人版放在星球中&#xff0c;我们开始正文。 大纲&#xff1a; 直接依赖接口依赖…

vue 腾讯云 javascript sdk + 简单富文本组件设计+实战

<template><div><quill-editor v-model"content" ref"myQuillEditor" :options"editorOption" change"onEditorChange"input"handleInput"></quill-editor><!-- 链接添加对话框 --><el-di…

【论文阅读笔记】In Search of an Understandable Consensus Algorithm (Extended Version)

1 介绍 分布式一致性共识算法指的是在分布式系统中&#xff0c;使得所有节点对同一份数据的认知能够达成共识的算法。且算法允许所有节点像一个整体一样工作&#xff0c;即使其中一些节点出现故障也能够继续工作。之前的大部分一致性算法实现都是基于Paxos&#xff0c;但Paxos…

前端Vue项目中腾讯地图SDK集成:经纬度与地址信息解析的实践

在前端开发中&#xff0c;我们经常需要将经纬度信息转化为具体的地址信息&#xff0c;这对于定位、地图展示等功能至关重要。Vue作为现代前端框架的代表&#xff0c;其组件化开发的特性使得我们能够更高效地实现这一功能。本文将介绍如何在Vue项目中集成腾讯地图SDK&#xff0c…

一个 基于nuxt3 + vite + ts 搭建的 网盘服务 (附带部署教程)

目录 介绍技术选型功能介绍代码地址部署安装 node 环境打包代码安装 pm2 去 后台运行代码安装一个nginx 介绍 最近 有个卖课的朋友 谈到 网盘没有目录分享的功能&#xff0c;我之前嫖了他太多课了&#xff0c;出于感激给他写个小服务。 在线地址&#xff1a; http://godboxs.c…

SpringMVC源码深度解析(上)

今天&#xff0c;聊聊SpringMVC框架的原理。SpringMVC属于Web框架&#xff0c;它不能单独存在&#xff0c;需要依赖Servlet容器&#xff0c;常用的Servlet容器有Tomcat、Jetty等&#xff0c;这里以Tomcat为例进行讲解。老规矩&#xff0c;先看看本项目的层级结构&#xff1a; 需…