fastspar微生物相关性推断

fastspar

简介

fastspar是基于Sparcc通过C编写的,速度更快,内存消耗更少。sparcc是基于OTU的原始count数,通过log转换和标准化去除传统相对丰度的天然负相关(因为所有OTU之和为1,某些OTU丰度高另外一些自然就少,导致最后出现正相关少负相关多的假象)。
FastSpar是SparCC算法的C++实现,比原来的Python2版本快几千倍,并且使用的内存少得多。FastSpar的实现提供了线程支持和一个P值估计器,该估计器考虑了重复数据排列的可能性(进一步的细节见本文)。

FastSpar目前正在开发中,可能缺乏完整程序中预期的某些功能。虽然FastSpar的工作一般没有问题,但必须特别注意确保提供正确格式的输入数据。

安装

Conda
To install through conda, use:

conda install -c bioconda -c conda-forge fastspar

使用

1.相关性推断

要运行FastSpar,你必须有BIOM tsv格式文件(没有元数据)中的绝对OTU计数。fake_data.tsv(来自最初的SparCC实现)将作为一个例子。

  • 输入文件:fake_data.tsv(fastspar/tests/data/fake_data.tsv,63k):制表符分隔,行为OTU,列为样本号
  • header(第一行第一列一定是#OTU ID,否则报错)
    在这里插入图片描述
    在这里插入图片描述
    计算相关性
(base) [yutao@myosin data]$ mkdir out; time fastspar --otu_table fake_data.tsv --correlation out/median_correlation.tsv --covariance out/median_covariance.tsv -t 30 &>1.log&
# 耗时2.4s
# -t <int>, --threads <int> Number of threads (default: 1)
-c <path>, --otu_table <path>OTU input OTU count table-r <path>, -correlation <path>Correlation output table-a <path>, --covariance <path>Covariance output table-i <int>, --iterations <int>Number of interations to perform (default: 50)
  • median_correlation.tsv
    在这里插入图片描述
    结果是一个相关系数对称矩阵,对角线是自身的相关系数为1,其他例如1行2列表示OTU0 vs OTU1相关系数为0.7265
  • median_covariance.tsv
    在这里插入图片描述
    迭代次数(SparCC相关性估计的轮次)和排除迭代次数(发现和排除高度相关的OTU对的次数)也可以进行调整。
fastspar --iterations 50 --exclude_iterations 20 --otu_table tests/data/fake_data.tsv --correlation median_correlation.tsv --covariance median_covariance.tsv

进一步,可以增加排除相关OTU对的最小阈值

fastspar --threshold 0.2 --otu_table tests/data/fake_data.tsv --correlation median_correlation.tsv --covariance median_covariance.tsv

2.精确P值的计算

有几种方法可以计算推断出的相关关系的p值。在这里,我们选择使用基于稳健互换的方法。这个过程包括从原始OTU计数数据的随机排列中推断出相关关系。每个p值的大小与随机排列的数据中观察到的更极端的相关性的频率有关。在下面的例子中,我们从1000个自举相关性中计算出p值。

1.首先我们生成1000个自举计数。

mkdir bootstrap_counts
(base) [yutao@myosin data]$ mkdir bootstrap_counts;  fastspar_bootstrap --otu_table tests/data/fake_data.tsv --number 1000 --prefix bootstrap_counts/fake_data
# 耗时1s
# --number 迭代次数
# -t 线程
  • 会在out下生成1000个重抽样列表
(base) [yutao@myosin data]$ ls out/
fake_data_0.tsv    fake_data_326.tsv  fake_data_552.tsv  fake_data_779.tsv
fake_data_100.tsv  fake_data_327.tsv  fake_data_553.tsv  fake_data_77.tsv
fake_data_101.tsv  fake_data_328.tsv  fake_data_554.tsv  fake_data_780.tsv
fake_data_102.tsv  fake_data_329.tsv  fake_data_555.tsv  fake_data_781.tsv

2.然后推断每个自举计数的相关性(在所有可用进程中并行运行)。
这里使用parallel进行并行计算

mkdir bootstrap_correlation
parallel fastspar --otu_table {} --correlation bootstrap_correlation/cor_{/} --covariance bootstrap_correlation/cov_{/} -i 5 ::: bootstrap_counts/*
# 1000次计算耗时14s
# bootstrap_correlation/cor_{/},bootstrap_correlation/cov_{/} 表示输出文件名是cor_,cov_加输入文件名
(base) [yutao@myosin data]$ ls bootstrap_correlation
cor_fake_data_0.tsv    cor_fake_data_701.tsv       cov_fake_data_400.tsv
cor_fake_data_100.tsv  cor_fake_data_702.tsv       cov_fake_data_401.tsv
cor_fake_data_101.tsv  cor_fake_data_703.tsv       cov_fake_data_402.tsv
cor_fake_data_102.tsv  cor_fake_data_704.tsv       cov_fake_data_403.tsv
cor_fake_data_103.tsv  cor_fake_data_705.tsv       cov_fake_data_404.tsv
cor_fake_data_104.tsv  cor_fake_data_706.tsv       cov_fake_data_405.tsv
cor_fake_data_105.tsv  cor_fake_data_707.tsv       cov_fake_data_406.tsv
  • 注意,此步骤需要有足够的存储,1278 * 85(大小400k)的矩阵产生了24G的中间结果

3.根据这些相关性,然后计算出p值

fastspar_pvalues --otu_table tests/data/fake_data.tsv --correlation median_correlation.tsv --prefix bootstrap_correlation/cor_fake_data_ --permutations 1000 --outfile pvalues.tsv
# 耗时2s
-c/--otu_table <path>OTU input table used to generated correlations-r/--correlation <path>Correlation table generated by FastSpar-p/--prefix <str>Prefix output of bootstrap output files-n/--permutations <int>Number of permutations/ bootstraps-o/--outfile <path>Output p-value matrix filename

线程
如果FastSpar是用OpenMP编译的,线程可以通过在命令行中调用–threads <thread_number>来使用几个工具。

fastspar --otu_table tests/data/fake_data.txt --correlation median_correlation.tsv --covariance median_covariance.tsv --iterations 50 --threads 10

解析成长表

# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
library(Hmisc) # 生成相关性/P-value矩阵,测试使用
library(dplyr) # 合并数据
setwd("C:/北大真菌细菌互作/互作")
tax <- read.table('Taxonomy_info.txt', header = T)
# 解析两个对角线矩阵,cormat是fastspac的相关性矩阵,pmat是fastspac的p-value矩阵
flattenCorrMatrix <- function(cormat, pmat){ut <- upper.tri(cormat) data.frame( row = rownames(cormat)[row(cormat)[ut]], column = rownames(cormat)[col(cormat)[ut]], cor =(cormat)[ut], p = pmat[ut] )
}r <- read.table('median_correlation.tsv', header = T, row.names = 1,check.names = F, comment.char = "", sep = "\t") #忽略注释字符
View(head(r))
pv <- read.table('pvalues.tsv', header = T, row.names = 1,check.names = F, comment.char = "", sep = "\t") #忽略注释字符
View(head(pv))res <- flattenCorrMatrix(r, pv)View(head(res))
write.table(res, 'Correlation_result.tsv', quote = F)dim(res)
filter <- subset(res, res$p < 0.05)
dim(filter)# 添加分组信息
View(head(tax))
r1 <- left_join(filter, tax, by = c('row' = 'ID'))
r1 <- left_join(r1, tax, by = c('column' = 'ID'))
# 添加互作类型
r1$Type <- 'null'r1$Type[r1$Kingdom.x == "Fungi" & r1$Kingdom.y == "Fungi"] <- "FF"
r1$Type[r1$Kingdom.x == "Bacteria" & r1$Kingdom.y == "Bacteria"] <- "BB"
r1$Type[r1$Kingdom.x == "Fungi" & r1$Kingdom.y == "Bacteria"] <- "FB"
r1$Type[r1$Kingdom.x == "Bacteria" & r1$Kingdom.y == "Fungi"] <- "FB"# 保留细菌vs真菌
write.table(r1, 'All_Correlation_result.tsv', quote = F, sep = "\t")data <- read.table('Bacteria_and_fungi_for_sparcc.tsv', header = T, row.names = 1, comment.char = "", sep = '\t' )
cor.test(as.matrix(data[1,]), as.matrix(data[2760,]), method = 'spearman')
#举个栗子
# res3 <- rcorr(as.matrix(mtcars[,1:7]))
# res <- flattenCorrMatrix(res3$r, res3$P)

报错

  • error while loading shared libraries: libmkl_rt.so
(base) [yutao@myosin data]$ fastspar_bootstrap --help
fastspar_bootstrap: error while loading shared libraries: libmkl_rt.so: cannot open shared object file: No such file or directory

解决

(base) [yutao@myosin data]$ conda install -c intel mkl
  • libc++abi: terminating with uncaught exception of type std::invalid_argument: stof: no conversion;Abort trap: 6
    要求输入的OTU表格为数值数据,除首行和首列外,其他均为数值,不能出现NA,对输入的大表数据可以先通过数据筛选确认符合上述情况。

Citation

fastspac github
If you use this tool, please cite the FastSpar paper and original SparCC paper:

Watts, S. C., Ritchie, S. C., Inouye, M., & Holt, K. E. (2018). FastSpar: rapid and scalable correlation estimation for compositional data. Bioinformatics. doi: 10.1093/bioinformatics/bty734

Friedman, J. & Alm, E.J. (2017). Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687.

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

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

相关文章

晶振分频【FPGA】

所有数据对齐晶振。 6分频&#xff1a;【1】 module divider_six // 6分频 【0~2】 ( input wire sys_clk , //系统时钟 50MHz input wire sys_rst_n , //全局复位 output reg clk_out //对系统时钟 6 分频后的信号 );reg [1:0] cnt; //用于计数的寄存器 //cnt:计数器从 0 到…

springboot定时服务

上一篇文章【修改定时时间&#xff0c;定时任务及时生效】 是定时任务与功能项目共用一个&#xff1b; 我目前所在公司的定时服务是专门有一个项目处理&#xff0c;然后定时查询库里面的定时信息配置。 话不多说&#xff0c;上程序 数据库设置 create table SCHEDULER_JOB…

一个“Hello, World”Flask应用程序

如果您访问Flask网站&#xff0c;会看到一个非常简单的示例应用程序&#xff0c;只有5行代码。为了不重复那个简单的示例&#xff0c;我将向您展示一个稍微复杂一些的示例&#xff0c;它将为您编写大型应用程序提供一个良好的基础结构。 应用程序将存在于包中。在Python中&…

AIGC:使用生成对抗网络GAN实现MINST手写数字图像生成

1 生成对抗网络 生成对抗网络&#xff08;Generative Adversarial Networks, GAN&#xff09;是一种非常经典的生成式模型&#xff0c;它受到双人零和博弈的启发&#xff0c;让两个神经网络在相互博弈中进行学习&#xff0c;开创了生成式模型的新范式。从 2017 年以后&#x…

论文浅尝 | 基于对多条思维链的元推理实现智能问答

笔记整理&#xff1a;屠铭尘&#xff0c;浙江大学硕士&#xff0c;研究方向为知识图谱 链接&#xff1a;https://arxiv.org/abs/2304.13007 1. 动机 1.1 Chain of Thought的诞生 尽管大语言模型在许多自然语言处理任务上表现出色&#xff0c;但由于其本质是token by token的类似…

使用easyui前端框架快速构建一个crud应用

本篇文章将会详细介绍jquery easyui前端框架的使用&#xff0c;通过创建一个crud应用来带大家快速掌握easyui的使用。 easyui是博主最喜欢的前端框架&#xff0c;没有之一&#xff0c;因为它提供了多种主题&#xff0c;而且有圆润的各种组件。 目录 一、快速开始 二、准备工作…

小程序如何设置下单提示语句

下单提示会展示在购物车和提交订单页面&#xff0c;它可以帮助商家告知客户事项&#xff0c;提高用户体验和减少错误操作。例如提示&#xff1a;商品是否包邮、某些区域是否发货、商品送达时间等等。 在小程序管理员后台->配送设置处&#xff0c;填写下单提示。在设置下单提…

Java自学第8课:电商项目(3) - 重新搭建环境

由于之前用的jdk和eclipse&#xff0c;以及mysql并不是视频教程所采用的&#xff0c;在后面运行源码和使用作者提供源码时&#xff0c;总是报错&#xff0c;怀疑&#xff1a; 1 数据库有问题 2 jdk和引入的jar包不匹配 3 其他什么未知的错误&#xff1f; 所以决定卸载jdk e…

SQL入门语句

MySQL和SQL的区别是什么&#xff1f;之间是什么关系&#xff1f; SQL&#xff08;Structured Query Language&#xff09;是用于管理和操作关系型数据库&#xff08;RDBMS&#xff09;的标准语言。SQL还可以用于这些RDBMS&#xff1a;MySQL、Oracle、Microsoft SQL Server、Pos…

Kafka中遇到的错误:

1、原因&#xff1a;kafka是一个去中心化结果的&#xff0c;所以在启动Kafka的时候&#xff0c;每一个节点上都需要启动。 启动的命令&#xff1a;kafka-server-start.sh -daemon /usr/local/soft/kafka_2.11-1.0.0/config/server.properties

AMESim 2021安装教程

主要是AMESim的安装 写在前面&#xff0c;由于项目需要&#xff0c;需要自学AMESim&#xff0c;因此需要安装这个软件&#xff0c;目前仅仅安装使用&#xff0c;还不涉及到与MATLAB的联合仿真&#xff0c;老板说用 RT LAB半实物仿真平台&#xff0c;但是简单搜了一下&#xff0…

Flutter笔记:绘图示例 - 一个简单的(Canvas )时钟应用

Flutter笔记 绘图示例 - 一个简单的&#xff08;Canvas &#xff09;时钟应用 作者&#xff1a;李俊才 &#xff08;jcLee95&#xff09;&#xff1a;https://blog.csdn.net/qq_28550263 邮箱 &#xff1a;291148484163.com 本文地址&#xff1a;https://blog.csdn.net/qq_2855…

如何记录血压的波动情况

import pandas as pd from plotnine import * import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [Microsoft YaHei] 记录时间(time)、收缩压(SBP)、舒张压(DBP)&#xff1a; df pd.DataFrame({ time: [2023-11-01 08:30, 2023-11-02 21:00, 2023-11-0…

前端-第一部分-HTML

一.初识HTML 1.1 HTML 简介 HTML 全称为 HyperText Mark-up Language&#xff0c;翻译为超文本标签语言&#xff0c;标签也称作标记或者元素。HTML 是目前网络上应用最为广泛的技术之一&#xff0c;也是构成网页文档的主要基石之一。HTML文本是由 HTML 标签组成的描述性文本&a…

力扣最热一百题——每日温度

Python后面的文章&#xff0c;内容都比较多&#xff0c;但是同时我又想保持每天更新的速度&#xff0c;所以Python的文章我继续打磨打磨&#xff0c;先更新一篇算法的文章。 一身正气报国家&#xff0c;旁无乱境不恋她 ヾ(◍∇◍)&#xff89;&#xff9e; 力扣题号&#xff1a…

css呼吸效果实现

实现一个图片有规律的大小变化&#xff0c;呈现呼吸效果&#xff0c;怎么用CSS实现这个呼吸效果呢 一.实现 CSS实现动态效果可以使用动画( animation)来属性实现&#xff0c;放大缩小效果可以用transform: scale来实现&#xff0c;在这基础上有了动画&#xff0c;就可以设置一个…

ps人像怎么做渐隐的效果?

photoshop怎么制作人像渐隐的图片效果&#xff1f;渐隐效果需要使用渐变来实现&#xff0c;下面我们就来看看详细的教程。 首先&#xff0c;我们打开Photoshop&#xff0c;点击屏幕框选的【打开】&#xff0c;打开一张背景图片。 下面&#xff0c;我们点击左上角【文件】——【…

Exploration by random network distillation论文笔记

Exploration by Random Network Distillation (2018) 随机网络蒸馏探索 0、问题 这篇文章提出的随机网络蒸馏方法与Curiosity-driven Exploration by Self-supervised Prediction中提出的好奇心机制的区别&#xff1f; 猜想&#xff1a;本文是基于随机网络蒸馏提出的intrin…

楼宇天台视频AI智能监管方案,时刻保障居民安全

一、背景需求分析 我们经常能看到这样的新闻报道&#xff0c;小孩登上小区的天台玩耍&#xff0c;因为家长和物业人员发现得晚&#xff0c;没有及时制止&#xff0c;结果导致意外事故的发生。此前&#xff0c;在某小区就有居民拍下多名儿童在小区高层住宅的楼顶玩耍跳跃&#…

Pytorch R-CNN目标检测-汽车car

概述 目标检测(Object Detection)就是一种基于目标几何和统计特征的图像分割,它将目标的分割和识别合二为一,通俗点说就是给定一张图片要精确的定位到物体所在位置,并完成对物体类别的识别。其准确性和实时性是整个系统的一项重要能力。 R-CNN的全称是Region-CNN(区域卷积神经…