数据分析:基于sparcc的co-occurrence网络

介绍

Sparcc是基于16s或metagenomics数据等计算组成数据之间关联关系的算法。通常使用count matrix数据。

安装Sparcc软件

git clone git@github.com:JCSzamosi/SparCC3.git
export PATH=/path/SparCC3:$PATHwhich SparCC.py

导入数据

注:使用rarefy抽平的count matrix数据

library(dplyr)
library(tibble)dat <- read.table("dat_rarefy10000_v2.tsv", header = T)

过滤数据

filter_fun <- function(prof=dat , tag="dat ", cutoff=0.005){# prof=dat # tag="dat " # cutoff=0.005dat <- cbind(prof[, 1, drop=F], prof[, -1] %>% summarise(across(everything(), function(x){x/sum(x)}))) %>%column_to_rownames("OTUID")#dat.cln <- dat[rowSums(dat) > cutoff, ]remain <- apply(dat, 1, function(x){length(x[x>cutoff])}) %>% data.frame() %>%setNames("Counts") %>%rownames_to_column("OTUID") %>%mutate(State=ifelse(Counts>1, "Remain", "Discard")) %>%filter(State == "Remain")# countcount <- prof %>% filter(OTUID%in%remain$OTUID)filename <- paste0("../dataset/Sparcc/", tag, "_rarefy10000_v2_", cutoff, ".tsv")write.table(count, file = filename, quote = F, sep = "\t", row.names = F)# relative abundancerelative <- dat %>% rownames_to_column("OTUID") %>%filter(OTUID%in%remain$OTUID)filename <- paste0("../dataset/Sparcc/", tag, "_rarefy10000_v2_", cutoff, "_rb.tsv")write.table(relative, file = filename, quote = F, sep = "\t", row.names = F)
}filter_fun(prof=dat, tag="dat", cutoff=0.005)
filter_fun(prof=dat, tag="dat", cutoff=0.001)

Result: 保留过滤后的count matrix和relative abundance matrix两类型矩阵

sparcc analysis

该过程分成两步:1.计算相关系数;2.permutation test计算p值.

  • iteration 参数使用default -i 20

  • permutation 参数: 1000次

# Step 1 - Compute correlations
python /data/share/toolkits/SparCC3/SparCC.py sxtr_rarefy10000_v2_0.001.tsv -i 20 --cor_file=sxtr_sparcc.tsv > sxtr_sparcc.log
echo "Step 1 - Compute correlations Ended successfully!"# Step 2 - Compute bootstraps
python /data/share/toolkits/SparCC3/MakeBootstraps.py sxtr_rarefy10000_v2_0.001.tsv -n 1000 -t bootstrap_#.txt -p pvals/ >> sxtr_sparcc.log
echo "Step 2 - Compute bootstraps Ended successfully!"# Step 3 - Compute p-values
for n in {0..999}; do /data/share/toolkits/SparCC3/SparCC.py pvals/bootstrap_${n}.txt -i 20 --cor_file=pvals/bootstrap_cor_${n}.txt >> sxtr_sparcc.log; done
python /data/share/toolkits/SparCC3/PseudoPvals.py sxtr_sparcc.tsv pvals/bootstrap_cor_#.txt 1000 -o pvals/pvals.two_sided.txt -t two_sided >> sxtr_sparcc.log
echo "Step 3 - Compute p-values Ended successfully!"# step 4 - Rename file
mv pvals/pvals.two_sided.txt sxtr_pvals.two_sided.tsv
mv cov_mat_SparCC.out sxtr_cov_mat_SparCC.tsv
echo "step 4 - Rename file Ended successfully!"

co-occurrence network

网络图要符合以下要求:

  1. 保留相互之间显著差异(p < 0.05)OTU;

  2. genus分类学水平表示OTU来源;

  3. OTU间相关性用颜色区分,且线条粗细代表相关系数大小;

  4. OTU点大小表示其丰度大小;

  5. 统计网络中正负相关数目;

导入画图数据

library(igraph)
library(psych)dat_cor <- read.table("dat_cov_mat_SparCC.tsv", header = T, row.names = 1)
dat_pval <- read.table("dat_pvals.two_sided.tsv", header = T, row.names = 1)
dat_rb5 <- read.table("dat_rarefy10000_v2_0.005_rb.tsv", header = T, row.names = 1)
dat_tax <- read.csv("dat_taxonomy.csv")

画图

  • 数据处理

  • 数据可视化

  • 数据存储

cornet_plot <- function(mcor=dat_cor, mpval=dat_pval, mrb=dat_rb5, tax=dat_tax, type="dat_5"){# mcor <- dat_cor# mpval <- dat_pval# mrb <- dat_rb5# tax <- dat_tax# type="dat_05"# trim all NA in pvalue < 0.05mpval[mpval > 0.05] <- NAremain <- apply(mpval, 1, function(x){length(x[!is.na(x)])}) %>% data.frame() %>%setNames("counts") %>%rownames_to_column("OTUID") %>%filter(counts > 0)remain_pval <- mpval[as.character(remain$OTUID), as.character(remain$OTUID)]# remove non significant edges remain_cor <- mcor[as.character(remain$OTUID), as.character(remain$OTUID)]for(i in 1:nrow(remain_pval)){for(j in 1:ncol(remain_pval)){if(is.na(remain_pval[i, j])){remain_cor[i, j] <- 0}}}# OTU relative abundance and taxonomy rb_tax <- mrb %>% rownames_to_column("OTUID") %>%filter(OTUID%in%as.character(remain$OTUID)) %>%group_by(OTUID) %>%rowwise() %>%mutate(SumAbundance=mean(c_across(everything()))) %>%ungroup() %>%inner_join(tax, by="OTUID") %>%dplyr::select("OTUID", "SumAbundance", "Genus") %>%mutate(Genus=gsub("g__Candidatus", "Ca.", Genus),Genus=gsub("_", " ", Genus)) %>%mutate(Genus=factor(as.character(Genus)))# 构建igraph对象igraph <- graph.adjacency(as.matrix(remain_cor), mode="undirected", weighted=TRUE, diag=FALSE)# 去掉孤立点bad.vs <- V(igraph)[degree(igraph) == 0]igraph <- delete.vertices(igraph, bad.vs)# 将igraph weight属性赋值到igraph.weightigraph.weight <- E(igraph)$weight# 做图前去掉igraph的weight权重,因为做图时某些layout会受到其影响E(igraph)$weight <- NAnumber_cor <- paste0("postive correlation=", sum(igraph.weight > 0), "\n","negative correlation=",  sum(igraph.weight < 0))# set edge color,postive correlation 设定为red, negative correlation设定为blueE.color <- igraph.weightE.color <- ifelse(E.color > 0, "red", ifelse(E.color < 0, "blue", "grey"))E(igraph)$color <- as.character(E.color)# 可以设定edge的宽 度set edge width,例如将相关系数与edge width关联E(igraph)$width <- abs(igraph.weight)# set vertices sizeigraph.size <- rb_tax %>% filter(OTUID%in%V(igraph)$name) igraph.size.new <- log((igraph.size$SumAbundance) * 1000000)V(igraph)$size <- igraph.size.new# set vertices colorigraph.col <- rb_tax %>% filter(OTUID%in%V(igraph)$name)pointcolor <- c("green","deeppink","deepskyblue","yellow","brown","pink","gray","cyan","peachpuff")pr <- levels(igraph.col$Genus)pr_color <- pointcolor[1:length(pr)]levels(igraph.col$Genus) <- pr_colorV(igraph)$color <- as.character(igraph.col$Genus)# 按模块着色# fc <- cluster_fast_greedy(igraph, weights=NULL)# modularity <- modularity(igraph, membership(fc))# comps <- membership(fc)# colbar <- rainbow(max(comps))# V(igraph)$color <- colbar[comps]filename <- paste0("../figure/03.Network/", type, "_Sparcc.pdf")pdf(file = filename, width = 10, height = 10)plot(igraph,main="Co-occurrence network",layout=layout_in_circle,edge.lty=1,edge.curved=TRUE,margin=c(0,0,0,0))legend(x=.8, y=-1, bty = "n",legend=pr,fill=pr_color, border=NA)text(x=.3, y=-1.2, labels=number_cor, cex = 1.5)dev.off()# calculate OTU remain_cor_sum <- apply(remain_cor, 1, function(x){res1 <- as.numeric(length(x[x>0]))res2 <- as.numeric(length(x[x<0]))res <- c(res1, res2)}) %>% t() %>% data.frame() %>%setNames(c("Negative", "Positive")) %>%rownames_to_column("OTUID")file_cor <- paste0("../figure/03.Network/", type, "_Sparcc_negpos.csv")write.csv(remain_cor_sum, file = file_cor, row.names = F)
}

运行画图函数

cornet_plot(mcor=dat_cor, mpval=dat_pval, mrb=dat_rb5, tax=dat_tax, type="dat_5")

在这里插入图片描述

R information

sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)Matrix products: defaultlocale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    
system code page: 936attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     other attached packages:
[1] psych_2.0.9  igraph_1.2.5 tibble_3.0.3 dplyr_1.0.2 loaded via a namespace (and not attached):[1] lattice_0.20-41  crayon_1.3.4     grid_4.0.2       R6_2.5.0         nlme_3.1-150     lifecycle_0.2.0  magrittr_1.5    [8] pillar_1.4.6     rlang_0.4.7      rstudioapi_0.11  vctrs_0.3.4      generics_0.1.0   ellipsis_0.3.1   tools_4.0.2     
[15] glue_1.4.2       purrr_0.3.4      parallel_4.0.2   xfun_0.19        yaml_2.2.1       compiler_4.0.2   pkgconfig_2.0.3 
[22] mnormt_2.0.2     tmvnsim_1.0-2    knitr_1.30       tidyselect_1.1.0

参考

  1. SparCC3

  2. sparcc.pdf

  3. sparcc tutorial

  4. Co-occurrence网络图在R中的实现

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

如若内容造成侵权/违法违规/事实不符,请联系瞎胡闹网进行投诉反馈,一经查实,立即删除!

相关文章

牛客小白月赛93

B交换数字 题目&#xff1a; 思路&#xff1a;我们可以知道&#xff0c;a*b% mod (a%mod) * (b%mod) 代码&#xff1a; void solve(){int n;cin >> n;string a, b;cin >> a >> b;for(int i 0;i < n;i )if(a[i] > b[i])swap(a[i], b[i]);int num1…

图片无损压缩工具-VIKY

一、前言 Viky v3.4是一款功能强大的图片压缩工具&#xff0c;它能够提供高效的图片无损压缩服务。通过使用独特的压缩算法&#xff0c;该软件在显著减小图片文件大小的同时&#xff0c;还保持了图像的清晰度和色彩饱和度&#xff0c;确保了图像质量的优异表现。 二、软件特点…

AS-VJ900实时视频拼接系统产品介绍:两画面视频拼接方法和操作

目录 一、实时视频拼接系统介绍 &#xff08;一&#xff09;实时视频拼接的定义 &#xff08;二&#xff09;无缝拼接 &#xff08;三&#xff09;AS-VJ900功能介绍 1、功能 2、拼接界面介绍 二、拼接前的准备 &#xff08;一&#xff09;摄像机选择 &#xff08;二&a…

区块链 | NFT 水印:Review on Watermarking Techniques(三)

&#x1f34d;原文&#xff1a;Review on Watermarking Techniques Aiming Authentication of Digital Image Artistic Works Minted as NFTs into Blockchains 一个 NFT 的水印认证协议 可以引入第三方实体来实现对交易的认证&#xff0c;即通过使用 R S A \mathsf{RSA} RSA…

力扣每日一题37:解数独

目录 题目 大致思路 方法一&#xff1a;回溯剪枝&#xff08;正常人能想出来的&#xff09; 方法二&#xff1a;位运算优化剪枝 需要使用的位运算技巧 代码 位运算怎么就优化了呢&#xff1f; 方法三&#xff1a;枚举优化。 官解代码 方法四&#xff1a;舞蹈链算法 题…

FreeRTOS任务调度器

目录 1、什么是任务调度器 2、FreeRTOS中的任务调度器 2.1 抢占式调度 2.2 时间片调度 2.3 协作式调度 3、任务调度案例分析 3.1 实验需求 3.2 CubeMX配置 3.3 代码实现 3.3.1 uart.c 重定向printf 3.3.2 打开freertos.c并添加代码 3.3.4 代码现象 1、什么是任务调度…

7 系列 FPGA 产品介绍及选型

目录 Spartan-7 FPGAsArtix-7 FPGAsKintex-7 FPGAsVirtex-7 FPGAsFPGA芯片命名规则DSP资源BRAM资源Transceivers 资源Transceivers 总带宽I/O 个数及带宽参考文档 Spartan-7 FPGAs Artix-7 FPGAs Kintex-7 FPGAs Virtex-7 FPGAs FPGA芯片命名规则 DSP资源 BRAM资源 Transceiver…

day5Qt作业

服务器端 #include "widget.h" #include "ui_widget.h"Widget::Widget(QWidget *parent): QWidget(parent), ui(new Ui::Widget) {ui->setupUi(this);//准备组件&#xff0c;初始化组件状态this->setFixedSize(800,600);chatwidget new QListWidge…

【Linux】Linux——Centos7安装Tomcat

1.下载Tomcat 安装包 官网地址&#xff1a;Apache Tomcat - Apache Tomcat 9 Software Downloadshttps://tomcat.apache.org/download-90.cgi 2.将下载的安装包上传到 Xftp 上&#xff0c;我是直接放到 usr 下了 3.将安装包解压到 /usr/local/ tar -zxvf apache-tomcat-9.0.8…

Java+SpringBoot+JSP实现在线心理评测与咨询系统

前言介绍 随着互联网技术的高速发展&#xff0c;人们生活的各方面都受到互联网技术的影响。现在人们可以通过互联网技术就能实现不出家门就可以通过网络进行系统管理&#xff0c;交易等&#xff0c;而且过程简单、快捷。同样的&#xff0c;在人们的工作生活中&#xff0c;也就…

249 基于matlab的MED、OMEDA、MOMEDA、MCKD信号处理方法

基于matlab的MED、OMEDA、MOMEDA、MCKD信号处理方法。最小熵反褶积(MED)&#xff0c;最优最小熵反卷积调整卷积 (OMEDA),多点最优最小熵解卷积调整&#xff08;Multipoint Optimal Minimum Entropy Deconvolution Adjusted&#xff0c;MOMEDA&#xff09;&#xff0c;最大相关峭…

Neuralink首个脑机接口患者:打游戏、搞研究两不误,重获自主能力

今年1月28日&#xff0c;Neuralink首次将侵入式脑机接口植入人类患者Noland Arbaugh的大脑。100天后&#xff0c;这家由埃隆马斯克创立的公司公布了最新的进展。从Neuralink的更新中我们可以看到&#xff0c;Arbaugh的恢复情况超出预期&#xff0c;他的用户体验也非常积极。 原…

Android 按钮Button点击音效

一、新建工程 编译运行&#xff0c;确保工程无误&#xff0c;这里不过多赘述。 二、UI布局 添加两个播放音效Button <?xml version"1.0" encoding"utf-8"?> <LinearLayout xmlns:android"http://schemas.android.com/apk/res/android"…

芋道系统springcloud模块启动报错,枚举类不能为空

问题描述&#xff1a; Error starting ApplicationContext. To display the conditions report re-run your application with debug enabled. 2024-05-10 15:50:15.756 | ERROR 9120 | main [TID: N/A] o.s.b.d.LoggingFailureAnalysisReporter | ************************…

腾讯云coding代码托管平台配置问题公钥拉取失败提示 Permission denied(publickey)

前言 最近在学校有个课设多人开发一个游戏&#xff0c;要团队协作&#xff0c;选用了腾讯云的coding作为代码管理仓库&#xff0c;但在配置的时候遇到了一些问题&#xff0c;相比于github&#xff0c;发现腾讯的coding更难用&#xff0c;&#xff0c;&#xff0c;这里记录一下…

【漫画版】指挥官的排序战术:快速排序算法解密

作者介绍&#xff1a;10年大厂数据\经营分析经验&#xff0c;现任字节跳动数据部门负责人。 会一些的技术&#xff1a;数据分析、算法、SQL、大数据相关、python&#xff0c;欢迎探讨交流 欢迎加入社区&#xff1a;码上找工作 作者专栏每日更新&#xff1a; LeetCode解锁1000题…

stm32单片机学习路线

第一步 编程及硬件基础知识 1.掌握C语言基础 作为STM32的主要编程语言&#xff0c;C语言的基础知识是必不可少的。建议通过书籍、在线课程或者教学视频系统地学习C语言的基础知识&#xff0c;包括语法、数据类型、函数、指针等。 2.了解电子电路基础 对于单片机开发来说&…

面向对象 03:类与对象的创建、初始化和使用,通过 new 关键字调用构造方法,以及创建对象过程的内存分析

一、前言 记录时间 [2024-05-10] 系列文章简摘&#xff1a; Java 笔记 01&#xff1a;Java 概述&#xff0c;MarkDown 常用语法整理 Java 笔记 11&#xff1a;Java 方法相关内容&#xff0c;方法的设计原则&#xff0c;以及方法的定义和调用 面向对象 01&#xff1a;Java 面向对…

ESP32引脚入门指南(四):从理论到实践(PWM)

引言 ESP32 作为物联网领域的明星微控制器&#xff0c;除了强大的Wi-Fi和蓝牙功能&#xff0c;还内置了丰富的外设资源&#xff0c;其中就包括高级的PWM&#xff08;脉冲宽度调制&#xff09;功能。本文将深入探讨ESP32的PWM引脚&#xff0c;解析其工作原理&#xff0c;并通过…

2024年数维杯B题完整代码和思路论文讲解与分析

2024数维杯数学建模完整代码和成品论文已更新&#xff0c;获取↓↓↓↓↓ https://www.yuque.com/u42168770/qv6z0d/bgic2nbxs2h41pvt?singleDoc# 2024数维杯数学建模B题45页论文和代码已完成&#xff0c;代码为全部问题的代码 论文包括摘要、问题重述、问题分析、模型假设、…