文章MSM_metagenomics(五):共现分析

欢迎大家关注全网生信学习者系列:

  • WX公zhong号:生信学习者
  • Xiao hong书:生信学习者
  • 知hu:生信学习者
  • CDSN:生信学习者2

介绍

本教程是使用一个Python脚本来分析多种微生物(即strains, species, genus等)的共现模式。

数据

大家通过以下链接下载数据:

  • 百度网盘链接:https://pan.baidu.com/s/1f1SyyvRfpNVO3sLYEblz1A
  • 提取码: 请关注WX公zhong号_生信学习者_后台发送 复现msm 获取提取码

Python packages required

  • pandas >= 1.3.5
  • matplotlib >= 3.5.0
  • seaborn >= 0.11.2

Co-presence pattern analysis

使用step_curve_drawer.py 做共线性分析

  • 代码
#!/usr/bin/env python"""
NAME: step_curve_drawer.py
DESCRIPTION: This script is to analyze the co-prsense of multiple species in different categories,by drawing step curves.
"""import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import sys
import argparse
import textwrapdef read_args(args):# This function is to parse argumentsparser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter,description = textwrap.dedent('''\This program is to do draw step curves to analyze co-presense of multiple species in different groups.'''),epilog = textwrap.dedent('''\examples:step_curve_drawer.py --abundance_table <abundance_table_w_md.tsv> --variable <variable_name> --species_number <nr_sps> --output <output.svg>'''))parser.add_argument('--abundance_table',nargs = '?',help = 'Input the MetaPhlAn4 abundance table which contains only a group of species one wants to analyze their co-presense state, with metadata being wedged.',type = str,default = None)parser.add_argument('--variable',nargs = '?',help = 'Specify the header of the variable in the metadata table you want to assess. For example, \[Diet] variable columns has three categries - [vegan]/[Flexitarian]/[Omnivore].',type = str,default = None)parser.add_argument('--minimum_abundance',nargs = '?',help = 'Specify the minimum abundance used for determining presense. note: [0, 100] and [0.0] by default',type = float,default = 0.0)parser.add_argument('--species_number',nargs = '?',help = 'Specify the total number of multiple species in the analysis.',type = int)parser.add_argument('--output',nargs = '?',help = 'Specify the output figure name.',type = str,default = None)parser.add_argument('--palette',nargs = '?',help = 'Input a tab-delimited mapping file where values are group names and keys are color codes.',type = str,default = None)return vars(parser.parse_args())class PandasDealer:"""This is an object for dealing pandas dataframe."""def __init__(self, df_):self.df_ = df_def read_csv(self):# Ths fucntion will read tab-delimitted file into a pandas dataframe.return pd.read_csv(self.df_, sep = '\t', index_col = False, low_memory=False)def rotate_df(self):# this function is to rotate the metaphlan-style table into tidy dataframe to ease searching work,df_ = self.read_csv()df_rows_lists = df_.values.tolist()rotated_df_dict = {df_.columns[0]: df_.columns[1:]}for i in df_rows_lists:rotated_df_dict[i[0]] = i[1:]rotated_df = pd.DataFrame.from_dict(rotated_df_dict)return rotated_dfclass CopEstimator:def __init__(self, sub_df_md):self.sub_df_md = sub_df_md # sub_df_md: a subset of dataframe which contains only a group of species one wants to do co-presence analysis.def make_copresense_df(self, factor, total_species_nr, threshold = 0.0):# factor: the factor you want to assess the category percentage.# total_species_nr: specify the total number of species you want to do co-presense analysis.rotated_df = PandasDealer(self.sub_df_md)rotated_df = rotated_df.rotate_df()cols = rotated_df.columns[-total_species_nr: ].to_list() categories = list(set(rotated_df[factor].to_list()))copresense = []cate_name = []ratios = []for c in categories:sub_df = rotated_df[rotated_df[factor] == c]species_group_df = sub_df[cols]species_group_df = species_group_df.apply(pd.to_numeric)species_group_df['total'] = species_group_df[cols].gt(threshold).sum(axis=1)for i in range(1, total_species_nr + 1):ratio = count_non_zero_rows(species_group_df, i)copresense.append(i)cate_name.append(c)ratios.append(ratio)return pd.DataFrame.from_dict({"copresense": copresense,factor: cate_name,"percentage": ratios})def count_non_zero_rows(df_, nr):total_rows = len(df_.index)sub_df = df_[df_['total'] >= nr]ratio = len(sub_df.index)/total_rowsreturn ratioclass VisualTools:def __init__(self, processed_df, factor):self.processed_df = processed_dfself.factor = factordef step_curves(self, opt_name, palette = None):categories = list(set(self.processed_df[self.factor].to_list()))if palette:palette_dict = {i.rstrip().split('\t')[0]: i.rstrip().split('\t')[1] for i in open(palette).readlines()}for c in categories:sub_df = self.processed_df[self.processed_df[self.factor] == c]plt.step(sub_df["percentage"]*100, sub_df["copresense"], label = c, color = palette_dict[c])else:for c in categories:sub_df = self.processed_df[self.processed_df[self.factor] == c]plt.step(sub_df["percentage"]*100, sub_df["copresense"], label = c)plt.title("Number of species in an individual if present")plt.xlabel("Percentage")plt.ylabel("Co-presense")plt.legend(title = self.factor)plt.savefig(opt_name, bbox_inches = "tight")if __name__ == "__main__":pars = read_args(sys.argv)cop_obj = CopEstimator(pars['abundance_table'])p_df = cop_obj.make_copresense_df(pars['variable'], pars['species_number'], pars['minimum_abundance'])vis_obj = VisualTools(p_df, pars['variable'])vis_obj.step_curves(pars['output'], palette = pars['palette'])
  • 用法
usage: step_curve_drawer.py [-h] [--abundance_table [ABUNDANCE_TABLE]] [--variable [VARIABLE]] [--minimum_abundance [MINIMUM_ABUNDANCE]] [--species_number [SPECIES_NUMBER]] [--output [OUTPUT]][--palette [PALETTE]]This program is to do draw step curves to analyze co-presense of multiple species in different groups.optional arguments:-h, --help            show this help message and exit--abundance_table [ABUNDANCE_TABLE]Input the MetaPhlAn4 abundance table which contains only a group of species one wants to analyze their co-presense state, with metadata being wedged.--variable [VARIABLE]Specify the header of the variable in the metadata table you want to assess. For example, [Diet] variable columns has three categries - [vegan]/[Flexitarian]/[Omnivore].--minimum_abundance [MINIMUM_ABUNDANCE]Specify the minimum abundance used for determining presense. note: [0, 100] and [0.0] by default--species_number [SPECIES_NUMBER]Specify the total number of multiple species in the analysis.--output [OUTPUT]     Specify the output figure name.--palette [PALETTE]   Input a tab-delimited mapping file where values are group names and keys are color codes.examples:python step_curve_drawer.py --abundance_table <abundance_table_w_md.tsv> --variable <variable_name> --species_number <nr_sps> --output <output.svg>

为了演示step_curve_drawer.py的使用,我们将绘制基于metaphlan相对丰度表特定于Segatalla copri(之前称为Prevotella copri)的八个谱系:./data/mpa4_pcopri_abundances_md.tsv的共现模式,这些数据来自MSMNon-MSM人群。MSMNon-MSM样本将使用自定义颜色进行标记,颜色分配来自一个颜色映射文件color map file: ./data/copresence_color_map.tsv

python step_curve_drawer.py \--abundance_table mpa_pcopri_abundances_md.tsv \--variable sexual_orientation \--species_number 8 \--palette copresence_color_map.tsv \--output copresence_plot.png

请添加图片描述

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

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

相关文章

持续集成jenkins+gitee

首先要完成gitee部署&#xff0c;详见自动化测试git的使用-CSDN博客 接下来讲如何从git上自动拉取代码&#xff0c;实现jenkins无人值守&#xff0c;定时执行测试&#xff0c;生成测试报告。 需要这三个安装包 由于目前的jenkins需要至少java11到java17的版本&#xff0c;所以…

JVM 性能分析案列——使用 JProfiler 工具分析 dump.hprof 堆内存快照文件排查内存溢出问题

在 windows 环境下实现。 参考文档 一、配置 JVM 参数 配置两个 JVM 参数&#xff1a; -XX:HeapDumpOnOutOfMemoryError&#xff0c;配置这个参数&#xff0c;会在发生内存溢出时 dump 生成内存快照文件&#xff08;xxx.hprof&#xff09;-XX:HeapDumpPathF:\logs&#xff…

力控算法每日一练:209. 长度最小的子数组(java)

给定一个含有 n 个正整数的数组和一个正整数 target 。 找出该数组中满足其总和大于等于 target 的长度最小的 子数组 [numsl, numsl1, ..., numsr-1, numsr] &#xff0c;并返回其长度。如果不存在符合条件的子数组&#xff0c;返回 0 。 class Solution {public int minSu…

代码随想录算法训练营第三十八天| 509. 斐波那契数 ,70. 爬楼梯,746. 使用最小花费爬楼梯

509. 斐波那契数 - 力扣&#xff08;LeetCode&#xff09; class Solution {public int fib(int n) {if (n < 1) {return n;}int[] dp new int[n 1];dp[0] 0;dp[1] 1;for (int i 2; i < n; i) {dp[i] dp[i - 1] dp[i - 2];}return dp[n];} } 70. 爬楼梯 - 力扣&am…

十二星座女、具有哪些情感特质。

白羊座&#xff08;奋不顾身&#xff09;。金牛座&#xff08;爱财如命&#xff09;。双子座&#xff08;灵活多变&#xff09;。 巨蟹座&#xff08;似水柔情&#xff09;。狮子座&#xff08;光明磊落&#xff09;。处女座&#xff08;尽善尽美&#xff09;。 天秤座&#xf…

安装wsl

安装wsl 先决条件&#xff1a; 打开控制面板->选择程序与功能->选择启动或关闭windows功能&#xff0c;将以下框选的勾选上 二、到Mircosoft store下载Ubuntu 三、如果以上都勾选了还报以下错误 注册表错误 0x8007019e Error code: Wsl/CallMsi/REGDB_E_CLASSNOTREG…

【three.js】旋转、缩放、平移几何体

目录 一、缩放 二、平移 三、旋转 四、居中 附源码 BufferGeometry通过.scale()、.translate()、.rotateX()、.rotateY()等方法可以对几何体本身进行缩放、平移、旋转,这些方法本质上都是改变几何体的顶点数据。 我们先创建一个平面物体,样子是这样的。 一、缩放 // 几何…

重新安装 Windows 10 后如何恢复丢失的数据?

“嗨&#xff0c;我的 Windows 10 崩溃了&#xff0c;所以我不得不重新安装它。我使用 USB 可启动驱动器重新安装了操作系统。但是&#xff0c;重新安装后&#xff0c;C 盘上的所有先前文件都丢失了。有什么方法可以恢复丢失的文件吗&#xff1f;” - Jacky 在大多数情况下&am…

如何在两个不同的conda环境中实现jupyter notebook共同使用,避免重复下载

前提&#xff1a;有2个conda环境&#xff0c;yes和py38_pytorch 其中&#xff0c;yes已经安装了jupyter notebook;py38_pytorch没有jupyter notebook 现在&#xff0c;实现在py38_pytorch用jupyter notebook 步骤&#xff1a; 1、激活py38_pytorch conda activate py38_p…

中小学电子教材下载办法(202406最简单的)

官方版本 现在能阅读电子教材的官方网站挺多的&#xff0c;例如 人民教育出版社-电子教材&#xff0c;还有 国家中小学智慧教育平台 &#xff0c;其他还有很多可在阅读的网站。由于平台的原因不能直接贴链接&#xff0c;大家可以通过搜索关键词找到网站。 如何下载 据我所知…

基于CentOS Stream 9平台安装MySQL8.4.0 LTS

1. 安装之前 1.1 查看系统版本 [rootcoisini /]# cat /etc/redhat-release CentOS Stream release 9 1.2 查看cpu架构 [rootcoisini /]# lscpu 架构&#xff1a; x86_64 CPU 运行模式&#xff1a; 32-bit, 64-bit 2. MySQL官方下载https://dev.mysql.com/downloads/mysql/ 或…

相亲交友APP系统|婚恋交友社交软件|语音聊天平台定制开发

在现代社会&#xff0c;婚恋交友已经成为了人们日常生活中的一项重要任务。为了方便用户进行相亲交友活动&#xff0c;各种相亲交友APP系统和婚恋交友社交软件应运而生。本文将介绍相亲交友APP系统、婚恋交友社交软件的开发以及语音聊天平台的定制开发的相关知识和指导。 一、…

special characters are not allowed

处理域名连接nacos读取配置异常 1 项目启动报错2 问题处理3 刷新依赖重启问题解决 1 项目启动报错 使用ip可以正在启动&#xff0c;但是使用域名报下面的错误 2024-06-15 17:37:22.981 ERROR 29268 --- [ main] c.a.c.n.c.NacosPropertySourceBuilder : parse …

餐厅点餐系统的设计

管理员账户功能包括&#xff1a;系统首页&#xff0c;个人中心&#xff0c;管理员管理&#xff0c;商品管理&#xff0c;用户管理&#xff0c;店家管理&#xff0c;广告管理 店家账户功能包括&#xff1a;系统首页&#xff0c;个人中心&#xff0c;商品管理&#xff0c;广告管…

牛客小白月赛96 解题报告 | 珂学家

前言 题解 A. 最少胜利题数 签到 n1 len(set(input())) n2 len(set(input()))if n1 < n2:n1, n2 n2, n1print (-1 if n1 6 else n1 - n2 1)B. 最少操作次数 思路: 分类讨论 只有-1,0,1,2这四种结果 特判 01, 10 n int(input()) s input()# 枚举 from collectio…

Windows10 MySQL(8.0.37)安装与配置

一、MySQL8.0.37下载 官网下载链接&#xff1a; https://dev.mysql.com/downloads/ 解压文件&#xff0c;解压到你想要的位置 二、新建MySQL配置文件 右键新建文本文档 新建my.txt文件 编辑my.txt文件&#xff0c;输入以下内容 [mysqld] # 设置 3306 端口 port3306 # 设…

SQLServer使用 PIVOT 和 UNPIVOT行列转换

在SQL Server中&#xff0c;PIVOT是一个用于将行数据转换为列数据的操作。它特别适用于将多个行中的值转换为多个列的情况&#xff0c;并在此过程中执行聚合操作。以下是关于SQL Server中PIVOT操作的详细解释和示例&#xff1a; 1、本文内容 概述语法备注关键点简单 PIVOT 示…

15.RedHat认证-Ansible自动化运维(上)

15.RedHat认证-Ansible自动化运维(上) RHCE8-RH294 Ansible自动化&#xff08;Ansible版本是2.8.2&#xff09; Ansible介绍 1.Ansible是什么&#xff1f; Ansible是一个简单的强大的无代理的自动化运维工具&#xff08;Ansible是自动化运维工具&#xff09;Ansible特点 简…

RPC知识

一、为什么要有RPC&#xff1a; HTTP协议的接口&#xff0c;在接口不多、系统与系统交互较少的情况下&#xff0c;解决信息孤岛初期常使用的一种通信手段&#xff1b;优点就是简单、直接、开发方便&#xff0c;利用现成的HTTP协议进行传输。 但是&#xff0c;如果是一个大型的网…

[大模型]XVERSE-7B-chat FastAPI 部署

XVERSE-7B-Chat为XVERSE-7B模型对齐后的版本。 XVERSE-7B 是由深圳元象科技自主研发的支持多语言的大语言模型&#xff08;Large Language Model&#xff09;&#xff0c;参数规模为 70 亿&#xff0c;主要特点如下&#xff1a; 模型结构&#xff1a;XVERSE-7B 使用主流 Deco…