Decontam去污染:一个尝试

为了程序运行的便利性,不想将Decontam放到windows的Rstudio里面运行,需要直接在Ubuntu中运行,并且为了在Decontam时进行其他操作,使用python去运行R

 首先你需要有一个conda环境,安装了R,Decontam,phyloseq,rpy2

mamba create -c bioconda -c conda-forge -p {os.path.join(args['conda_env_dir'], 'LBSWrap-Decontam')} r-base=4.3.1 python=3.10 -y >/dev/null 2>&1
mamba install -p {os.path.join(args['conda_env_dir'], 'LBSWrap-Decontam')} bioconda::bioconductor-decontam -y >/dev/null 2>&1
mamba install -p {os.path.join(args['conda_env_dir'], 'LBSWrap-Decontam')} bioconda::bioconductor-phyloseq -y >/dev/null 2>&1
mamba install -p {os.path.join(args['conda_env_dir'], 'LBSWrap-Decontam')} biopython -y >/dev/null 2>&1
os.path.join(args['conda_env_dir'], 'LBSWrap-Decontam', 'pip')} install rpy2 >/dev/null 2>&1

然后建立这个python文件,chmod 755给权限就行了

需要有你的fasta文件(用于输出通过过滤的和未通过过滤的fasta序列),otu文件,metadata文件 

#! /usr/bin/env python
#########################################################
# run Decontam (method="prevalence") in Ubuntu through python
# written by PeiZhong in IFR of CAAS

import argparse
import rpy2.robjects as robjects
from Bio import SeqIO
import os

parser = argparse.ArgumentParser(description='run Decontam (method="prevalence") in Ubuntu through python')
parser.add_argument('--otu_txt', '-otu', type=str, required=True, help='< OTU table file >')
parser.add_argument('--metadata_txt', '-metadata', required=True, type=str, help='< Metadata file >')
parser.add_argument('--way', '-w', type=str, required=True, help='< choice from isNotContaminant or isContaminant >')
parser.add_argument('--threshold', '-t', type=float, default=None, help='< Threshold , isNotContaminant default = 0.5, isContaminant default = 0.1 >')
parser.add_argument('--result_txt', '-r', type=str, required=True, help='< result txt from Decontam >')
parser.add_argument('--input_fasta', '-ifa', type=str, required=True, help='< your fasta file >')
parser.add_argument('--output_fasta', '-ofa', type=str, required=True, help='< clean fasta according to the result of Decontam >')
parser.add_argument('--contamination_fasta', '-cfa', type=str, required=True, help='< contamination fasta according to the result of Decontam >')

args = parser.parse_args()

for filepath in [args.otu_txt, args.metadata_txt, args.input_fasta]:
    if not os.path.exists(filepath):
        raise FileNotFoundError(f"The file {filepath} does not exist.")

if args.way not in ["isNotContaminant", "isContaminant"]:
    raise ValueError("The 'way' parameter must be either 'isNotContaminant' or 'isContaminant'.")

if args.threshold is None:
    args.threshold = 0.5 if args.way == "isNotContaminant" else 0.1

shi, fou = ("TRUE", "FALSE") if args.way == "isNotContaminant" else ("FALSE", "TRUE")

r_code = f"""
library(ggplot2)
library(decontam)
library(phyloseq)

otu <- read.table("{args.otu_txt}", header=TRUE, row.names = 1)
sample <- read.table("{args.metadata_txt}", header=TRUE, row.names = 1)

otu_table <- otu_table(t(otu), taxa_are_rows = FALSE)
sample_data <- sample_data(sample)
ps <- phyloseq(otu_table, sample_data)

sample_data(ps)$is.neg <- sample_data(ps)$Sample_or_Control == "Control"
contamdf.notcontam <- {args.way}(ps, method="prevalence", neg="is.neg", threshold={args.threshold},detail=TRUE) # isNotContaminant -> True means non-contaminants
write.table(contamdf.notcontam, file="{args.result_txt}", sep="\t", quote=F, row.names=T)
"""
robjects.r(r_code)

df, df_contam = {}, {}
with open(args.result_txt,"r") as read_result:
    for line in read_result.readlines():
        if "p.prev" not in line and len(line.split("\t")) >= 6:
            name, determine = line.split("\t")[0].strip("\n"), line.split("\t")[6].strip("\n")
            if determine == shi:
                df[name] = 0
            elif determine == fou:
                df_contam[name] = 0

with open(args.input_fasta, 'r') as fasta_file:
    sequences = SeqIO.to_dict(SeqIO.parse(fasta_file, 'fasta'))

clean_sequences = {key: sequences[key] for key in df}
contaminated_sequences = {key: sequences[key] for key in df_contam}

with open(args.output_fasta, 'w') as clean_file, open(args.contamination_fasta, 'w') as contaminated_file:
    SeqIO.write(clean_sequences.values(), clean_file, 'fasta')
    SeqIO.write(contaminated_sequences.values(), contaminated_file, 'fasta')

OTU示例和metadata示例参考....emmm算了比较麻烦,我给你们看我的文件

注意R里面"-"这个横杠会被识别为".",容易出错

使用示例

/home/zhongpei/hard_disk_sda2/zhongpei/Software/my_script/Decontam.py -t 0.5 \ 
-otu ${zubie}_${ruanjian}_decontam_pre.txt -metadata ${zubie}_metadata.txt \
 -w isContaminant -r ${zubie}_${ruanjian}_decontam.report \ 
-ifa ${zubie}_${ruanjian}_all.fa -ofa ${zubie}_${ruanjian}_decontam.fa \ 
-cfa ${zubie}_${ruanjian}_notgood.fa

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

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

day 49 动态规划 part 10● 121. 买卖股票的最佳时机 ● 122.买卖股票的最佳时机II

看了题解&#xff0c;第一种暴力&#xff0c;两个for循环。 class Solution { public:int maxProfit(vector<int>& prices) {int result 0;for (int i 0; i < prices.size(); i) {for (int j i 1; j < prices.size(); j){result max(result, prices[j] -…

使用scrapy爬取蜻蜓FM

创建框架和项目 ### 1. 创建虚拟环境 conda create -n spiderScrapy python3.9 ### 2. 安装scrapy pip install scrapy2.8.0 -i https://pypi.tuna.tsinghua.edu.cn/simple### 3. 生成一个框架并进入框架 scrapy startproject my_spider cd my_spider### 4. 生成项目 scrapy …

LeetCode:143.重排链表

143. 重排链表 解题过程 /*** Definition for singly-linked list.* public class ListNode {* int val;* ListNode next;* ListNode() {}* ListNode(int val) { this.val val; }* ListNode(int val, ListNode next) { this.val val; this.next next; …

数据结构——堆的应用 Topk问题

&#x1f49e;&#x1f49e; 前言 hello hello~ &#xff0c;这里是大耳朵土土垚~&#x1f496;&#x1f496; &#xff0c;欢迎大家点赞&#x1f973;&#x1f973;关注&#x1f4a5;&#x1f4a5;收藏&#x1f339;&#x1f339;&#x1f339; &#x1f4a5;个人主页&#x…

实验一:华为VRP系统的基本操作

1.1实验介绍 1.1.1关于本实验 本实验通过配置华为设备&#xff0c;了解并熟悉华为VRP系统的基本操作 1.1.2实验目的 理解命令行视图的含义以及进入离开命令行视图的方法 掌握一些常见的命令 掌握命令行在线帮助的方法 掌握如何撤销命令 掌握如何使用命令快捷键 1.1.3实验组网 …

挑战杯 基于设深度学习的人脸性别年龄识别系统

文章目录 0 前言1 课题描述2 实现效果3 算法实现原理3.1 数据集3.2 深度学习识别算法3.3 特征提取主干网络3.4 总体实现流程 4 具体实现4.1 预训练数据格式4.2 部分实现代码 5 最后 0 前言 &#x1f525; 优质竞赛项目系列&#xff0c;今天要分享的是 基于深度学习机器视觉的…

React组件(函数式组件,类式组件)

函数式组件 <!DOCTYPE html> <html lang"en"> <head><meta charset"UTF-8"><meta name"viewport" content"widthdevice-width, initial-scale1.0"><title>React Demo</title> <!-- 引…

阿里云服务器9元1个月优惠价格表

阿里云服务器9元1个月优惠价格表&#xff0c;用不上9元&#xff0c;又降价了&#xff0c;只要5元。阿里云服务器一个月多少钱&#xff1f;最便宜5元1个月。阿里云轻量应用服务器2核2G3M配置61元一年&#xff0c;折合5元一个月&#xff0c;2核4G服务器30元3个月&#xff0c;2核2…

深入探讨AI团队的角色分工

目录 前言1. 软件工程师&#xff1a;构建系统基石的关键执行者2. 机器学习工程师&#xff1a;数据与模型的塑造专家3. 机器学习研究员&#xff1a;引领算法创新的智囊4. 机器学习应用科学家&#xff1a;理论与实践的巧妙连接5. 数据分析师&#xff1a;洞察数据&#xff0c;智慧…

C语言初学10:typedef

一、作用 为用户定义的数据类型取一个新名字 二、对结构体使用typedef定义新的数据类型名字 #include <stdio.h> #include <string.h>typedef struct Books //使用 typedef 来定义一个新的数据类型名字 {char title[50];} book;int main( ) {//book是typedef定…

图片在div完全显示

效果图&#xff1a; html代码&#xff1a; <div class"container" style" display: flex;width: 550px;height: 180px;"><div class"box" style" color: red; background-color:blue; width: 50%;"></div><div …

基于冠豪猪优化算法(Crested Porcupine Optimizer,CPO)的无人机三维路径规划(MATLAB)

一、无人机路径规划模型介绍 无人机三维路径规划是指在三维空间中为无人机规划一条合理的飞行路径&#xff0c;使其能够安全、高效地完成任务。路径规划是无人机自主飞行的关键技术之一&#xff0c;它可以通过算法和模型来确定无人机的航迹&#xff0c;以避开障碍物、优化飞行…

企业如何安全参与开源项目?

【开源三句半】 企业参与开源潮&#xff0c; 安全创新都重要&#xff0c; 持续投入不可少&#xff0c; 眼光独到。 开源已经成为构建现代软件的常见方式&#xff0c;这不仅局限于IT技术本身&#xff0c;更推动了多个行业的数字化发展。企业决定引入开源项目打造商业软件时&…

NeRF模型NeRF模型

参考视频&#xff1a;https://www.youtube.com/watch?vHfJpQCBTqZs&ab_channelVision%26GraphicsSeminaratMIT NeRF模型的输入输出: 输入: (x, y, z): 一个三维空间坐标,代表场景中的一个位置点(θ, φ): 视线方向,θ表示与y轴的夹角,φ表示与x轴的夹角,用两个角度可以…

random标准模块

一、概述 在 Python 中&#xff0c;random 是一个内置模块&#xff0c;用于生成随机数。它提供了各种用于生成随机数的函数&#xff0c;包括伪随机数生成器、随机序列操作等。 1、需要导包 不会自动导入&#xff0c;需要显示的将random模块导入 import random2、源码分析&…

力扣最热100题——56.合并区间

吾日三省吾身 还记得梦想吗 正在努力实现它吗 可以坚持下去吗 目录 吾日三省吾身 力扣题号&#xff1a;56. 合并区间 - 力扣&#xff08;LeetCode&#xff09; 题目描述 Java解法一&#xff1a;排序然后原地操作 具体代码如下 Java解法二&#xff1a;new一个list&#xf…

K8S - 在任意node里执行kubectl 命令

当我们初步安装玩k8s &#xff08;master 带 2 nodes&#xff09; 时 正常来讲kubectl 只能在master node 里运行 当我们尝试在某个 node 节点来执行时&#xff0c; 通常会遇到下面错误 看起来像是访问某个服务器的8080 端口失败了。 原因 原因很简单 , 因为k8s的各个组建&…

使用Tokeniser估算GPT和LLM服务的查询成本

将LLM集成到项目所花费的成本主要是我们通过API获取LLM返回结果的成本&#xff0c;而这些成本通常是根据处理的令牌数量计算的。我们如何预估我们的令牌数量呢&#xff1f;Tokeniser包可以有效地计算文本输入中的令牌来估算这些成本。本文将介绍如何使用Tokeniser有效地预测和管…

IM6ULL学习总结(四-七-1)输入系统应用编程

第7章 输入系统应用编程 7.1 什么是输入系统 ⚫ 先来了解什么是输入设备&#xff1f; 常见的输入设备有键盘、鼠标、遥控杆、书写板、触摸屏等等,用户通过这些输入设备与 Linux 系统进行数据交换。 ⚫ 什么是输入系统&#xff1f; 输入设备种类繁多&#xff0c;能否统一它们的…

解决 cx-programmer 梯形图中繁体中文乱码问题

我的情况 cx-programmer9.5是繁体版&#xff0c;梯形图编辑区中打出的字体&#xff0c;简体繁体 都是乱码。 但是状态栏显示注解是正常的繁体。 原因 简体和繁体的编码不一样。繁体的BIG5和简体的GB2312不能互转&#xff0c;A编码的用B解码也是乱码。 解决 把系统字体调整为繁…
最新文章