学习一下ARG-OAP的python代码

江有枫xx

已于 2023-11-26 08:34:40 修改

阅读量288

点赞数

文章标签：学习

于 2023-09-14 22:40:57 首次发布

本文链接：https://blog.csdn.net/youfeng_xjy/article/details/132880345

版权

Step1. 认识该软件使用的db

1)src/args_oap/db/ko30.fasta

19,951个氨基酸序列，都是些核糖体蛋白

2) src/args_oap/db/ko30_structure.txt

共19,950行，给ko30.fasta的序列名分配了一个BA号

3）gg85 （green gene 85%）

5,088条DNA序列

Step2. 看看stage_one

1）def init（运行前的基本检查）

比如有没有input文件夹，output文件夹里原来有没有东西，检查数据库文件是否齐全...看看就好。

这一部分在def __init__中

2）功能函数

def count_16s(self, file)

使用bwa（预过滤）和blastn（后过滤）提取16S（GreenGenes 16S rRNA数据库85%）读数。

i) 使用bwa预过滤

subprocess.run(['bwa', 'mem', '-t', str(self.thread), '-o', _tmp_16s_sam, settings._gg85, file], check=True, stderr=subprocess.DEVNULL)

self.thread：bwa的线程数
settings._gg85GreenGenes 16S rRNA 数据库的路径
check=True：这个参数表示在命令执行过程中，如果返回的退出状态码不为零，会引发一个CalledProcessError异常。换句话说，如果命令执行出错，将会引发异常。
stderr=subprocess.DEVNULL：将标准错误输出重定向到subprocess.DEVNULL，即将其丢弃，不保存到任何地方。

iI) 将 SAM 格式转换为 FASTA 格式

with open(_tmp_16s_fa, 'w') as f:
subprocess.run(['samtools', 'fasta', '-F4', '-F0x900', _tmp_16s_sam], check=True, stderr=subprocess.DEVNULL, stdout=f)

将结果写入指定的文件 _tmp_16s_fa 中。samtools fasta 命令用于此转换过程。选项 -F4 和 -F0x900 用于过滤未匹配的 reads 和二级质量不合格的 reads。

iii）使用 blastn 进行后续筛选

mt_mode = '1' if simple_count(_tmp_16s_fa)[0] / self.thread >= 2500000 else '0'
subprocess.run([
'blastn',
'-db', settings._gg85,
'-query', _tmp_16s_fa,
'-out', _tmp_16s_txt,
'-outfmt', ' '.join(['6']+settings.cols),
'-evalue', str(self.e1),
'-max_hsps', '1',
'-max_target_seqs', '1',
'-mt_mode', mt_mode,
'-num_threads', str(self.thread)], check=True, stderr=subprocess.DEVNULL)

这段代码使用 blastn 命令对转换后的 FASTA 文件进行进一步筛选。blastn 命令用于将查询序列与数据库中的序列进行比对。具体选项和参数包括：

-db：指定要比对的数据库，这里是 settings._gg85。
-query：指定要进行比对的查询序列文件，这里是 _tmp_16s_fa。
-out：指定输出比对结果的文件路径，这里是 _tmp_16s_txt。
-outfmt：指定输出结果的格式，这里使用的是 Blast tabular 格式，列名为 settings.cols 中的值。
-evalue：指定比对的 E-value 阈值。
-max_hsps 和 -max_target_seqs：分别指定每个查询序列的最大比对结果数和每个数据库序列的最大比对结果数。
-mt_mode 和 -num_threads：用于设定多线程模式和线程数。
mt_mode = '1' if simple_count(_tmp_16s_fa)[0] / self.thread >= 2500000 else '0' 也就是选择是设定多线程模式，关于mt_mode的解释如下
-mt_mode <Integer, (>=0 and =<1)>
Multi-thread mode to use in BLAST search:
0 (auto) split by database
1 split by queries
Default = `0'

iv) 处理比对结果并计算覆盖率

df = pd.read_table(_tmp_16s_txt, header=None, names=settings.cols)
if len(df)==0:
logger.warning("No 16S-like sequences found in file <{}>.".format(file))
return 0
else:
df['scov'] = df['length'] / df['slen']
if df['qseqid'].duplicated().sum()>0:
logger.warning('Duplicated qseqid in 16S.')
df = df[~df['qseqid'].duplicated()]
return df['scov'].sum()