学习一下ARG-OAP的python代码

Step1. 认识该软件使用的db

一共有这些个数据库

1)src/args_oap/db/ko30.fasta

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

出自 30 sets of ESCMGs，用30个基因（属于基本单拷贝基因家族）来估计微生物群落中基于宏基因组测序数据的细胞的平均基因组大小。一个群落的基因组大小将与该群落中必需单拷贝基因的相对丰度R成反比。

   引用文献自：Average genome size estimation improves comparative metagenomics and sheds light on the functional ecology of the human microbiome - PMC《平均基因组大小估计改进了比较宏基因组学，并揭示了人类微生物组的功能生态学》，这其实就是MicrobeCensus的文章

2) src/args_oap/db/ko30_structure.txt

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

 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()

def count_cells(self, file)

使用diamond计数基本单拷贝标记基因（细胞数）

df.groupby('ko30').apply(lambda x: sum(x['length'] / x['slen'])).sum() / 30

• df.groupby('ko30'):

  • 这部分代码将 df 按照 ko30 列的值进行分组。每个组将包含具有相同 ko30 值的所有行。

• .apply(lambda x: sum(x['length'] / x['slen'])):

  • 对每个分组，使用 apply 方法和一个匿名函数（lambda 函数）进行操作。
  • x 代表当前组的数据框（DataFrame）。
  • x['length'] / x['slen'] 计算每一行的 length 列与 slen 列的比值。这个比值可能表示某种比例或覆盖度。
  • sum(x['length'] / x['slen']) 计算当前组中所有这些比值的总和。

• .sum():

  • 这部分代码将所有组的结果相加，得到所有 groups 的比值总和。

• / 30:

  • 最后，将总和除以 30。这可能是因为在计算细胞数量时，假设每个 ko30 类别代表一个特定的标记基因，而你有 30 个这样的基因，因此用总和除以 30 来获取每个基因的平均值或比例。

def extract_seqs(self, file)

def run(self)

 3）运行

def run_stage_one(options):
    StageOne(vars(options)).run()

Step3. 看看stage_two

1)函数

def merge_files

功能：将提取的目标序列与元数据和结构信息合并，并根据不同的层级（如类型、亚型、基因）进行汇总

logger.info打印log文件

筛选：

self.qcov = self.qcov / 100 / 3 if self.dbtype == 'prot' else self.qcov / 100  # for prot need to lower qcov as qlen is in bp
self.length = self.length if self.dbtype == 'prot' else self.length * 3  # for nucl need to convert aa cut to bp


df['qcov'] = df['length'] / df['qlen']
df = df[(df['pident'] >= self.id) & (df['evalue'] <= self.e) & (df['length'] >= self.length) & (df['qcov'] >= self.qcov)]

nbps, nlines = simple_count(self.setting.extracted)：调用了一个名为simple_count的函数，并将返回的结果赋值给nbps和nlines两个变量。它用于计算文件中的碱基数和行数。

mt_mode = '1' if nbps / self.thread >= 2500000 else '0'：这是一个条件语句，根据nbps / self.thread的值是否大于等于2500000，将mt_mode变量设置为'1'或'0'。它用于决定一种处理模式。

根据self.dbtype的值是否为'prot'，将blast_mode变量设置为'blastx'或'blastn'

def extract_seqs

将提取的目标序列与metadata数据和structure文件合并，并根据不同的层级（type/subtype/gene）进行聚合。

df['qcov'] = df['length'] / df['qlen']

def run

如果self.blastout为空，那么会执行self.extract_seqs()函数来提取序列

Step4. 看看settings.py

定义一下类（比如Setting类，File类）啊方法啊什么的，比如指定self.setting.gg85的文件路径啥的

@property 是一个装饰器，用于将一个方法转换为属性

@dataclass 是 Python 的一个装饰器，自动为类添加一些特殊方法，比如 __init__、__repr__ 和 __eq__，使得类的实例化和比较更加方便。