单细胞|GeneTrajectory·基因轨迹

跑完了,记录一下,顺便写写我在使用中遇到的问题,欢迎讨论~

声明:我是用自己数据跑的,因为还未发表所以就还是借用官网的图啦~

1.准备

library(GeneTrajectory)
library(Seurat)
library(dplyr)
library(reticulate)
library(RColorBrewer)
library(ggplot2)
library(scales)

2.加载数据集

data_S <- readRDS("./yourscdata.rds")
DimPlot(data_S, group.by = "celltype", shuffle = T)

3.基因间距离计算 

选择基因

这里作者选择了前 2000 个高变基因中 1% 到 50% 的细胞表达的基因来演示基因间距离计算。具体可以看自己的选择,也可以选择一些你认为重要的基因。

assay <- "RNA"
DefaultAssay(data_S) <- assay
data_S <- FindVariableFeatures(data_S, nfeatures = 2000)
all_genes <- data_S@assays[[assay]]@var.features
expr_percent <- apply(as.matrix(data_S[[assay]]@data[all_genes, ]) > 0, 1, sum)/ncol(data_S)
genes <- all_genes[which(expr_percent > 0.01 & expr_percent < 0.5)]
length(genes)
## [1] 251

准备用于基因-基因距离计算的输入

# Compute the Diffusion Map cell embedding
data_S <- GeneTrajectory::RunDM(data_S)
# Calculate cell-cell graph distances over a cell-cell kNN graph
cell.graph.dist <- GetGraphDistance(data_S, K = 10)
# Coarse-grain the cell graph by grouping cells into `N`=500 "meta-cells"
cg_output <- CoarseGrain(data_S, cell.graph.dist, genes, N = 500)

计算基因间距离

使用 reticulate R 包设置 virtualenv,py_install安装gene-trajectory。

if(!reticulate::virtualenv_exists('gene_trajectory')){reticulate::virtualenv_create('gene_trajectory', packages=c('gene_trajectory'))
}
reticulate::use_virtualenv('gene_trajectory')
reticulate::py_install("gene-trajectory")
# Import the function to compute gene-gene distances
cal_ot_mat_from_numpy <- reticulate::import('gene_trajectory.compute_gene_distance_cmd')$cal_ot_mat_from_numpy
# Compute gene-gene distances 
gene.dist.mat <- cal_ot_mat_from_numpy(ot_cost = cg_output[["graph.dist"]], gene_expr = cg_output[["gene.expression"]], num_iter_max = 50000, show_progress_bar = TRUE)
rownames(gene.dist.mat) <- cg_output[["features"]]
colnames(gene.dist.mat) <- cg_output[["features"]]
dim(gene.dist.mat)

4.基因轨迹推断和可视化 

gene_embedding <- GetGeneEmbedding(gene.dist.mat, K = 5)$diffu.emb
gene_trajectory <- ExtractGeneTrajectory(gene_embedding, gene.dist.mat, N = 3, t.list = c(4,7,7), K = 5)
table(gene_trajectory$selected)
## 
## Trajectory-1 Trajectory-2 Trajectory-3 
##           43           55          153
# Visualize gene trajectories
par(mar = c(1.5,1.5,1.5,1.5))
scatter3D(gene_embedding[,1],gene_embedding[,2],gene_embedding[,3],bty = "b2", colvar = as.integer(as.factor(gene_trajectory$selected))-1,main = "trajectory", pch = 19, cex = 1, theta = 45, phi = 0,col = ramp.col(c(hue_pal()(3))))

可以得到基因的轨迹图,从得到的数据来看,每个gene都有一个order值。

5.可视化基因箱图 

 这一步比较耗时,自定义N.bin值,可以分开看每个轨迹的基因分布在细胞上的进展路径

# Seurat v4安装旧版本SeuratWrappers
# remotes::install_github('satijalab/seurat-wrappers@community-vignette')
library(SeuratWrappers)
data_S <- RunALRA(data_S)
data_S <- AddGeneBinScore(data_S, gene_trajectory, N.bin = 5, trajectories = 1:3, assay = "alra", reverse = c(F, F, T))# Visualize gene bin plots for each gene trajectory
FeaturePlot(data_S, pt.size = 0.05, features = paste0("Trajectory",1,"_genes", 1:5), ncol = 5, order = T) &scale_color_gradientn(colors = rev(brewer_pal(palette = "RdYlBu")(10))) & NoLegend() & NoAxes() & 
theme(title = element_text(size = 10))
FeaturePlot(data_S, pt.size = 0.05, features = paste0("Trajectory",2,"_genes", 1:5), ncol = 5, order = T) &scale_color_gradientn(colors = rev(brewer_pal(palette = "RdYlBu")(10))) & NoLegend() & NoAxes() & 
theme(title = element_text(size = 10))
FeaturePlot(data_S, pt.size = 0.05, features = paste0("Trajectory",3,"_genes", 1:5), ncol = 5, order = T) &scale_color_gradientn(colors = rev(brewer_pal(palette = "RdYlBu")(10))) & NoLegend() & NoAxes() & 
theme(title = element_text(size = 10))

如何比较不同条件的基因轨迹?

参考作者回复:

参考:Editing GeneTrajectory/scripts/mouse_dermal_example.R at main · KlugerLab/GeneTrajectory (github.com)

python版本可参考:KlugerLab/GeneTrajectory-python: Python implementation of Gene Trajectory (github.com) 

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

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

相关文章

OpenCV 入门(五) —— 人脸识别模型训练与 Windows 下的人脸识别

OpenCV 入门系列&#xff1a; OpenCV 入门&#xff08;一&#xff09;—— OpenCV 基础 OpenCV 入门&#xff08;二&#xff09;—— 车牌定位 OpenCV 入门&#xff08;三&#xff09;—— 车牌筛选 OpenCV 入门&#xff08;四&#xff09;—— 车牌号识别 OpenCV 入门&#xf…

STM32程序进入hardfault_handler()

背景&#xff1a; 假期前一直在修改代码&#xff0c;没有边改边测。节后回来测试代码&#xff0c;发现程序上电后很快就进入hardfault_handler&#xff08;&#xff09;中断。 导致程序反复复位。 查找原因&#xff1a; 在程序的_it.c文件里有几句代码&#xff0c;如果注释…

【陀螺仪JY61P维特智能】通过单片机修改波特率和角度参考的方法

根据官方文档&#xff1a; 修改波特率 1.解锁:FF AA 69 88 B5 1.1延时200ms 2.修改波特率:FF AA 04 06 00 2.1切换已修改的波特率然后重新发送解锁和保存指令 2.2解锁:FF AA 69 88 B5 2.3延时200ms 4.保存: FF AA 00 00 00 XY轴角度参考 角度参考是以传感器当前的实际位置&…

【系统分析师】系统分析部分

文章目录 1、系统分析概述2、详细调查2.1 为什么要做详细调查&#xff1f;2.2 详细调查的原则2.3 详细调查的内容2.4 详细调查的方法 3、现有系统分析3.1 获得系统的物理模型3.2 抽象出现有系统的逻辑模型3.3 建立新系统的逻辑模型3.4 建立新系统的物理模型 4、组织结构分析4.1…

记录汇川:电磁阀封装

二位电磁阀封装&#xff1a; 中封三位电磁阀封装&#xff1a; HMI&#xff1a;

用Redis延时队列搞定订单超时业务

Redis延时队列是一种用于在特定时间后执行任务的消息队列。它在许多场景中非常有用&#xff0c;比如订单超时自动关闭、定时提醒等。在Redis中&#xff0c;通常使用Sorted Set&#xff08;有序集合&#xff09;来实现延时队列&#xff0c;因为Sorted Set可以按照分数进行排序&a…

10页面结构分析

我们打开一个网页&#xff0c;都会有一个清晰的结构和布局上图中的标签就是用来划分各个部分区域用的。其中比较常用重要的是header、footer和nav&#xff0c;需要重点掌握。 下面是部分代码及效果演示 <header> <h2>网页头部</h2> </header><sec…

vue2结合element-ui实现TreeSelect 树选择功能

需求背景 在日常开发中&#xff0c;我们会遇见很多不同的业务需求。如果让你用element-ui实现一个 tree-select 组件&#xff0c;你会怎么做&#xff1f; 这个组件在 element-plus 中是有这个组件存在的&#xff0c;但是在 element-ui 中是没有的。 可能你会直接使用 elemen…

详解MySQL常用的数据类型

前言 MySQL是一个流行的关系型数据库管理系统&#xff0c;它支持多种数据类型&#xff0c;以满足不同数据处理和存储的需求。理解并正确使用这些数据类型对于提高数据库性能、确保数据完整性和准确性至关重要。本文将详细介绍MySQL中的数据类型&#xff0c;包括数值类型、字符…

大数据高级阶段面试题(实时)

1.Kafka的producer如何实现幂等性? ①开启幂等性&#xff0c;将Idempotent设置为true ②将ack设置为-1&#xff0c;确保相同的消息只会发送一次&#xff0c;避免重新发送 2.Kafka的ISR和OSR的作⽤分别是什么? ISR是副本和领导者的数据和状态要保持一致&#xff0c;如果出现…

游泳耳机哪个牌子好性价比高?这四款热榜游泳耳机必须要看!

随着生活品质的提升和健康意识的增强&#xff0c;游泳已成为许多人日常锻炼的首选。在享受水中畅游的同时&#xff0c;音乐成为了许多游泳爱好者的最佳伴侣。游泳耳机&#xff0c;作为一种专为水下运动设计的音频设备&#xff0c;近年来逐渐受到市场的青睐。然而&#xff0c;面…

新华三VRRP配置

新华三VRRP配置 配置步骤 (1).基础配置&#xff1a; CORE1&#xff1a; [CORE1]vlan 10 //创建vlan10 [CORE1-vlan10]int vlan 10 //进入vlanif 10 [CORE1-Vlan-interface10]ip add 192.168.10.1 24 //配置ip [CORE1-Vlan-interface10]int g1/0/2 //进入接口 [C…

【OCPP】ocpp1.6协议第3.11章节Reservations和第3.12章节Vendor-specific data transfer-介绍及翻译

目录 3.11章节Reservations 概述 3.11章节Reservations 译文 3.12章节Vendor-specific data transfer 概述 3.12章节Vendor-specific data transfer 译文 3.11章节Reservations 概述 OCPP1.6协议中的3.11章节关于“Reservations”主要对充电桩预定过程进行了定义和规定。 基…

一文了解CRM系统帮助中心:从认识到搭建

客户关系管理&#xff08;CRM&#xff09;系统是企业的一个重要部分。而CRM系统帮助中心为用户提供了便捷的支持服务&#xff0c;提升了用户体验&#xff0c;减少了企业运营成本。本文将从认识到搭建&#xff0c;带你全面了解CRM系统帮助中心。 一、认识CRM系统帮助中心 CRM系统…

品鉴中的艺术表达:如何将红酒与绘画、雕塑等艺术形式相结合

品鉴雷盛红酒不仅是一种味觉的享受&#xff0c;更是一种艺术的体验。将雷盛红酒与绘画、雕塑等艺术形式相结合&#xff0c;能够创造出与众不同的审美体验&#xff0c;进一步丰富品鉴的内涵。 首先&#xff0c;绘画作为视觉艺术的一种表现形式&#xff0c;能够通过色彩和构图来传…

Python爬虫实战:爬取小红书去水印图片

1. 思路分析 首先&#xff0c;在小红书 APP 中点击分享&#xff0c;获取到它的链接分享&#xff0c;如&#xff1a;www.xiaohongshu.com/discovery/i… 然后把它在浏览器中&#xff08;我用的是 chrome 浏览器&#xff09;打开。 按 F12 或者 Ctrl shift i 打开 开发者工具…

SG-8018CE晶体振荡器可编程规格书

SG-8018CE系列晶体振荡器是一个高性能、多功能且具有高度集成性的解决方案&#xff0c;它满足了现代电子系统的严格要求。其广泛的频率范围0.67 MHz到170 MHz&#xff0c;且频率调节精度达到1ppm&#xff0c;1.62 V至3.63V的宽广电源电压&#xff0c;使能&#xff08;OE&#x…

RTSP/Onvif安防监控系统EasyNVR级联视频上云系统EasyNVS报错“Login error”的原因排查与解决

EasyNVR安防视频云平台是旭帆科技TSINGSEE青犀旗下支持RTSP/Onvif协议接入的安防监控流媒体视频云平台。平台具备视频实时监控直播、云端录像、云存储、录像检索与回看、告警等视频能力&#xff0c;能对接入的视频流进行处理与多端分发&#xff0c;包括RTSP、RTMP、HTTP-FLV、W…

小程序预览或上传代码时,遇到app.json未找到某个wxml文件的解决方法

uniapp小程序&#xff0c;点击预览或者是上传代码&#xff0c;遇到app.json无法找到某个wxml文件的解决方法&#xff1a;清缓存 问题&#xff1a; message&#xff1a;Error: app.json: 未找到 ["subPackages"][3]["pages"][3] 对应的 subPackages4/pages/…

Rust 解决循环引用

导航 循环引用一、现象二、解决 循环引用 循环引用出现的一个场景就是你指向我&#xff0c;我指向你&#xff0c;导致程序崩溃 解决方式可以通过弱指针&#xff0c;而Rust中的弱指针就是Weak 在Rc中&#xff0c;可以实现&#xff0c;对一个变量&#xff0c;持有多个不可变引…