社区所有版块导航
Python
python开源   Django   Python   DjangoApp   pycharm  
DATA
docker   Elasticsearch  
aigc
aigc   chatgpt  
WEB开发
linux   MongoDB   Redis   DATABASE   NGINX   其他Web框架   web工具   zookeeper   tornado   NoSql   Bootstrap   js   peewee   Git   bottle   IE   MQ   Jquery  
机器学习
机器学习算法  
Python88.com
反馈   公告   社区推广  
产品
短视频  
印度
印度  
Py学习  »  Python

Ehrapy-基于python的临床多组学数据端到端分析工具(四)

灵活胖子的科研进步之路 • 1 周前 • 19 次点击  
教程首页

教程地址:https://ehrapy.readthedocs.io/en/stable/tutorials/notebooks/mimic_2_introduction.html

聚类分析

为了更好地理解数据,通常通过莱顿算法中实现的社区检测来确定集群是有用的。此外,聚类允许对集群之间发生变化的特征进行无偏检测,因此对我们来说很有用的。

聚类识别

ehrapy中的实现允许设置确定找到的集群数量的分辨率。调整参数通常很有用。

# 使用Leiden算法进行聚类分析
# 参数说明:
# - adata: AnnData对象,包含需要聚类的数据
# - resolution: 聚类分辨率参数,值为0.3
#   - 较小的resolution值会产生较少的聚类
#   - 较大的resolution值会产生较多的聚类
#   - 0.3是一个中等的分辨率,可以得到适中数量的聚类
# - key_added: 聚类结果存储的键名为"leiden_0_3"
#   - 聚类结果将保存在adata.obs["leiden_0_3"]中
#   - 每个细胞会被分配一个聚类标签
# 该聚类的作用:
# 1. 识别数据中的自然分组
# 2. 发现具有相似特征的患者群体
# 3. 为后续分析提供群体分类依据
ep.tl.leiden(adata, resolution=0.3, key_added="leiden_0_3")

leiden算法向存储集群的obs(leiden_0_3)添加了一个密钥。这些随后可以在UMAP嵌入中可视化。

adata.obs.head(4)
# 使用UMAP可视化Leiden聚类结果
# 参数说明:
# - adata: AnnData对象,包含UMAP降维和Leiden聚类结果
# - color: 指定用于着色的变量
#   - "leiden_0_3": 使用resolution=0.3的Leiden聚类标签进行着色
#   - 不同颜色代表不同的聚类群体
# - title: 设置图表标题为"Leiden 0.3"
#   - 表明这是resolution=0.3的Leiden聚类结果
# - size: 设置散点大小为20
# 该可视化的作用:
# 1. 直观展示Leiden聚类的结果
# 2. 观察不同聚类之间的分布和关系
# 3. 评估聚类效果的合理性
ep.pl.umap(adata, color=["leiden_0_3"], title="Leiden 0.3" , size=20)

接下来,我们可以探索某些特定于集群的特征,因此可以用于注释。

集群特性

为了识别特定于集群的标记,ehrapy提供了ep.tl.rank_features_groups()函数,该函数允许集群组之间的统计测试以确定显着丰富或降低的值

# 对Leiden聚类的不同组进行特征排序分析
# 参数说明:
# - adata: AnnData对象,包含数据和聚类结果
# - groupby: 指定用于分组的变量名
#   - "leiden_0_3": 使用resolution=0.3的Leiden聚类标签作为分组依据
#   - 会对每个聚类组进行特征排序
# 该分析的作用:
# 1. 识别每个聚类组的特征性标记
# 2. 找出区分不同组别的关键特征
# 3. 帮助理解每个聚类组的生物学意义
# 4. 为后续分析提供重要特征的参考
ep.tl.rank_features_groups(adata, groupby="leiden_0_3")
# 设置图形参数
# - figsize=(4, 4): 设置图形大小为4x4英寸
# - dpi=100: 设置图形分辨率为100dpi(每英寸点数)
ep.settings.set_figure_params(figsize=(4, 4), dpi=100)

# 可视化特征排序结果
# 参数说明:
# - adata: AnnData对象,包含特征排序的结果
# - key: 指定要可视化的特征排序结果的key
#   - "rank_features_groups": 使用之前计算的特征排序结果
# - ncols=2: 设置图形列数为2,即每行显示2个子图
# 该可视化的作用:
# 1. 展示每个聚类组中最具特征性的基因/特征
# 2. 通过热图形式直观显示特征表达模式
# 3. 帮助识别不同组别的标志性特征
# 4. 为生物学解释提供依据
ep.pl.rank_features_groups(adata, key="rank_features_groups", ncols=2)

我们还可以将每个集群的顶级功能作为DataFrame

# 获取特征排序的数据框
# - ep.ad.get_rank_features_df(): 从AnnData对象中提取特征排序结果
# - adata: 包含特征排序结果的AnnData对象
# - group=["0", "1", "2", "3", "4", "5"]: 指定要分析的聚类组别
df = ep.ad.get_rank_features_df(adata, group=["0""1""2""3""4""5"])

# 筛选显著性差异的特征
# - df.loc[]: 基于条件筛选数据框
# - df["logfoldchanges"] > 0: 选择表达量上调的特征(log fold change > 0)
# - df["pvals_adj"] 
# - &: 组合两个条件,同时满足才会被保留
df = df.loc[(df["logfoldchanges"] > 0) & (df["pvals_adj"
# 筛选第2组(group="2")的特征数据
# - df.loc[]: 基于条件定位数据框的行
# - df["group"] == "2": 选择group列值等于"2"的行
# - 返回一个包含第2组所有特征的数据框,包含:
#   - names: 特征名称
#   - scores: 特征得分
#   - logfoldchanges: 表达量变化倍数(log2)
#   - pvals: 原始p值
#   - pvals_adj: 校正后的p值
df.loc[df["group"] == "2",]

从这张表中,我们还可以提取每个集群中的顶部特征,并在UMAP上或作为小提琴图逐集群突出显示这些特征。

# 获取每个组别的前5个重要特征
# - df.groupby("group"): 按group列分组
# - .head(5): 每组取前5行数据
# - 这样可以获得每个聚类组最显著的5个特征
top_features = df.groupby("group").head(5)

# 提取唯一的特征名称并转换为Series
# - top_features["names"]: 获取names列
# - .unique(): 获取唯一值,去除重复
# - pd.Series(): 转换为pandas Series对象
# 这样可以得到所有组别中最重要特征的去重列表
top_features = pd.Series(top_features["names"].unique())

# 显示结果
# 输出所有组别中最重要的特征名称列表
top_features
# 设置图形参数
# - figsize=(3.8, 2): 设置图形大小为3.8x2英寸
# - dpi=100: 设置图形分辨率为100dpi
ep.settings.set_figure_params(figsize=(3.8, 2), dpi=100)

# 绘制第一组小提琴图
# - keys=["censor_flg", "mort_day_censored"]: 展示审查标志和死亡天数
# - groupby="leiden_0_3": 按leiden聚类结果分组
ep.pl.violin(adata, keys=["censor_flg""mort_day_censored"], groupby="leiden_0_3")

# 绘制第二组小提琴图
# - keys=["platelet_first", "age"]: 展示血小板首次检测值和年龄
# - groupby="leiden_0_3": 按leiden聚类结果分组
ep.pl.violin(adata, keys=["platelet_first""age"], groupby="leiden_0_3")

# 绘制第三组小提琴图
# - keys=["sapsi_first", "copd_flg"]: 展示SAPS评分和COPD标志
# - groupby="leiden_0_3": 按leiden聚类结果分组
ep.pl.violin(adata, keys=["sapsi_first""copd_flg"], groupby="leiden_0_3")

# 绘制第四组小提琴图
# - keys=["sofa_first", "liver_flg"]: 展示SOFA评分和肝功能标志
# - groupby="leiden_0_3": 按leiden聚类结果分组
ep.pl.violin(adata, keys=["sofa_first""liver_flg"], groupby="leiden_0_3")

集群注释

通过对集群特征的了解,以及上面的UMAP图,我们可以对集群进行注释。

# 初始化注释列,默认值为"NA"
adata.obs["annotation"] = "NA"

# 定义聚类组的注释字典
# - key: 聚类组编号(0-5)
# - value: 该组的特征描述
# - liver+/sofa+表示肝功能和SOFA评分较高
# - weight+表示体重较高
# - age+表示年龄较大
# - stroke+表示有中风
# - deceased+表示死亡率高
# - platelet+表示血小板计数高
# - malignancy+表示恶性肿瘤
# - copd+表示慢性阻塞性肺疾病
annotation = {
    "0""liver+/sofa+",
    "1""weight+"
    "2""age+/stroke+/deceased+",
    "3""platelet+",
    "4""age+/malignancy+/copd+/deceased+",
    "5""age+",
}

# 根据leiden聚类结果为每个样本添加注释
# - 遍历leiden_0_3列的值
# - 如果值存在于annotation字典中,使用对应的注释
# - 否则保持原始的聚类编号
adata.obs["annotation"] = [
    annotation[l] if l in annotation.keys() else l for l in adata.obs["leiden_0_3"]
]

附加下游分析

在这些基本的ehrapy分析步骤之后,可以执行额外的下游分析(另请参阅其他教程)。

PAGA

推断轨迹以了解动态过程和阶段转换也可能很感兴趣。ehrapy为此提供了几种轨迹推断算法。其中之一是基于分区的图抽象(PAGA)。

# 使用PAGA(Partition-based Graph Abstraction)算法分析聚类结果
# - adata: AnnData对象,包含数据和注释信息
# - groups="leiden_0_3": 使用leiden_0_3列作为分组依据
# PAGA算法可以:
# 1. 分析不同聚类之间的连接关系
# 2. 构建聚类的抽象图表示
# 3. 帮助理解数据的拓扑结构
ep.tl.paga(adata, groups="leiden_0_3")
# 绘制PAGA(Partition-based Graph Abstraction)可视化图
# - adata: AnnData对象,包含数据和注释信息
# - color: 指定用于着色的列名列表
#   - leiden_0_3: 使用Leiden聚类结果进行着色
#   - day_28_flg: 使用28天内死亡标志进行着色
# - cmap: 设置颜色映射为灰色到红色渐变
#   - 使用预定义的grey_red配色方案
# - title: 为两个子图设置标题
#   - "Leiden 0.3": 显示Leiden聚类结果
#   - "Died in less than 28 days": 显示28天内死亡情况
ep.pl.paga(
    adata,
    color=["leiden_0_3""day_28_flg"],
    cmap=ep.pl.Colormaps.grey_red.value,
    title=["Leiden 0.3""Died in less than 28 days"],
)
# 使用PAGA初始化位置计算UMAP降维
# - adata: AnnData对象,包含数据和注释信息
# - init_pos="paga": 使用PAGA的结果作为UMAP的初始位置
#   这样可以保持PAGA分析得到的全局拓扑结构
ep.tl.umap(adata, init_pos="paga")

# 绘制UMAP可视化图
# - adata: AnnData对象,包含数据和注释信息 
# - color=["annotation"]: 使用annotation列进行着色
#   annotation列包含了不同的分组标签
#   可以直观展示不同组别在UMAP空间中的分布
ep.pl.umap(adata, color=["annotation"])
# 使用force-directed graph layout算法绘制图形
# - adata: AnnData对象,包含数据和注释信息
# - init_pos="paga": 使用PAGA的结果作为初始位置
#   这样可以保持PAGA分析得到的全局拓扑结构
ep.tl.draw_graph(adata, init_pos="paga")

# 绘制force-directed graph可视化图
# - adata: AnnData对象,包含数据和注释信息
# - color: 指定用于着色的列名列表
#   - leiden_0_3: 使用Leiden聚类结果进行着色
#   - day_28_flg: 使用28天内死亡标志进行着色
# - legend_loc="on data": 将图例放置在数据点上
#   可以直观地展示不同组别的分布情况
ep.pl.draw_graph(adata, color=["leiden_0_3""day_28_flg"], legend_loc="on data")

导出结果

我们将所有计算和最终状态保存到h5ad文件中。然后可以使用ep.io.read()函数再次读取它,例如:

# 将AnnData对象保存为h5ad文件
# - "mimic_2.h5ad": 指定保存的文件名
#   - 使用h5ad格式,这是AnnData对象的标准存储格式
#   - 可以保存所有数据、注释和分析结果
# - adata: 要保存的AnnData对象
#   - 包含了原始数据
#   - 包含了所有的分析结果(如聚类、降维等)
#   - 包含了注释信息
# 保存后的文件可以通过ep.io.read()重新读取
ep.io.write("mimic_2.h5ad", adata)

结论

MIMIC-II IAC数据集包含来自1776个人的46个特征的电子健康记录(EHR)。这种高维数据不容易解释,当只关注选定的明确定义的特征时,可以监督许多有趣和以前未知的特征。为了克服这个障碍,我们在MIMIC-II IAC数据集上应用了ehrapy。

EHRapy基于AnnData数据结构和稀疏的管道,以实现高效分析。我们使用内置函数对数据进行预处理,通过缺失数据的插补执行QC并降低维度,从而产生PCA和UMAP降维分析。在执行所有这些步骤之后,我们通过可视化UMAP降维上的多个特征来探索数据,从而第一眼看到患者结构。为了以无偏见的方式识别患者群体,我们使用Leiden算法对数据进行聚类,从而产生7个不同的患者聚类。聚类特定特征的计算使我们能够根据最突出的标记对聚类进行注释。我们看到已故、年龄较高且患有严重合并症(如中风和COPD)的患者(聚类2+3)与具有较温和特征(如血小板和体重增加(聚类0+1)的患者之间存在强烈差异。靠近这两个集群的是另外两个集群,它们具有更严重的特征,如心率加快(集群5)和肝病SOFA评分高(集群6),表明潜在的患者轨迹。集群4与所有其他集群不同,由离开ICU几个月/几年后死亡的患者组成。

ep.print_versions()

-----
ehrapy              0.9.0
rich                NA
scanpy              1.10.4
session_info        1.0.0
-----
PIL                 11.0.0
anndata             0.11.1
array_api_compat    1.9.1
asttokens           NA
attr                24.3.0
autograd            NA
autograd_gamma      NA
cachetools          5.5.0
causallearn         NA
certifi             2024.12.14
charset_normalizer  3.4.0
cloudpickle         3.1.0
colorama            0.4.6
comm                0.2.2
cycler              0.12.1
cython_runtime      NA
dateutil            2.9.0.post0
db_dtypes           1.3.1
debugpy             1.8.11
...
Python 3.10.15 | packaged by Anaconda, Inc. | (main, Oct  3 2024, 07:22:19) [MSC v.1929 64 bit (AMD64)]
Windows-10-10.0.26100-SP0
-----
Session information updated at 2024-12-25 19:49
Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

最新课程-基于R语言的动态预测模型课程-胖子老师独自授课

开课目的及前言

预测模型作为真实世界研究的重要组成部分,其研究被广泛开展。但是,传统的预测模型利用基线数据对最终的生存结果进行预测,这种模型无法纳入患者在后续随访中可能会动态变化的重要数据(比如肿瘤标记物的动态变化)。 以上情况在统计学中会产生估计偏差情况,也是不符合临床实际的。近年来发展起来的动态预测模型方法,利用患者的多次随访数据,结合患者的基线数据,对最终患者的额生存结果(或类似的time to event事件)进行估计。其发文量呈现快速增长趋势。

在临床实际中,医生会根据患者的动态变化指标做出进一步诊断及治疗的判断。动态预测模型结合患者的纵向数据与最终的生存结果,对于最终结果进行更加准备的预测。由于当前R语言在医学统计工作中占据重要地位,但很多临床大夫、护士因为时间工作关系很难将R语言与临床科研相结合,故开设R语言动态预测模型课程,旨在快速让学员掌握统计工作中常用到的R语言,助力临床科研工作。天企助力(天津)生产力促进有限公司特举办“基于R语言的动态预测模型课程培训班”。

预测模型类文章目前总结起来发展经历了以下三个阶段:

  1. 基于传统流行病学的列线图模型(本质都是cox回归及glm回归),简单的统计学分析模型,是模型依赖的方法,临床上实际情况很难满足其前提假设,实际效果不好。

  2. 基于机器学习/深度学习的预测模型的构建(在数据上提高了维度,在算法上引入了机器学习),虽然算法上引入了机器学习模型,处理数据更加灵活,模型的假设也更少。但是在使用的数据上还是患者的一次基线数据进行预测,与临床实际不符。

  3. 基于纵向数据的动态预测模型(基于纵向多次随访数据,模型应用联合模型等动态预测模型方法),应用患者的多次随访数据对最终的生存结果进行预测,从数据和方法上都更类似于临床实际。

考虑到动态预测模型有以下特点,因此必然是后续高分文章的必备方法:

  1. 数据上必须有同一个患者的多次随访数据,相对于既往横断面一次基线数据,数据的收集难度更大,而且动态预测模型需拟合纵向的线性混合模型,因此需要的数据量较大。这就提示我们如果能收集到如上数据更加容易发高分文章。

  2. 应用方法学动态预测模型需首先掌握普通生存分析及普通预测模型的方法,并且还需要熟悉纵向数据分析的广义线性混合模型,再次基础上还需要掌握tidyverse语法基础来将自己的数据转换为满足函数要求的纵向数据,另外对于联合模型,模型的结合形式及变量选择也均需要从临床背景及统计学方法考虑。

近期高分文章举例

文章示例-动态预测模型预测筛查肠癌患者
文章示例-动态预测模型预测前列腺癌预后
文章示例-动态预测用于创伤外科
文章示例-动态预测对比传统模型在糖尿病患者中的应用
顶刊文章示例-动态预测模型用于肾移植后再次肾功能不全诊断
杂志情况

授课老师

灵活胖子-独自

双一流学校肿瘤学博士毕业,目前就职于国内五大肿瘤中心之一。科研方向为真实世界研究,生物信息学分析及人工智能研究。目前以第一或共同第一作者身份发表SCI论文10余篇,累计IF50+。目前与国内多个院校及医院有科研合作。联合翻译小组同学,在国内第一次将jmbayes2及dynamicLM全文翻译为中文并在公众号发表。

课程目录及安排

授课形式及时间

授课形式:远程在线实时直播授课。

授课时间:2024年12月开课,总课时不少于30小时,每周进行3-5小时的授课,有充分时间学习,预计6-8周完成所有授课内容。

答疑支持:建立课程专属微信群,1年内课程内容免费答疑。

视频回看:3年内免费无限次回看。

课程售价及售后保证

课程售价:总价3000元,报名可先交300元预定即可,开课后2周内交齐即可

对公转账等手续务必提前联系助教

承办公司:天企助力(天津)生产力促进有限公司

奖励政策:学员应用所学内容发表IF 10+文章可退还学费(具体要求及流程需要咨询助教)

报名咨询

可联系我的助教进行咨询

我的助教微信

助教联系电话:18502623993

正式通知

pdf版通知可联系助教获取


Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/177855
 
19 次点击