学习如何使用 Bash and/or Biopython 管理表格数据（来自 BLAST+6 格式）的资源

Question

我是运行 BLAST 并且想使用 BLAST+6 格式操作输出。

例如，我想获取每次命中的 E 值、查询覆盖率和标识，然后将它们代入一个方程，将这三者权衡成一个 "score"。然后我想获取所有分数并将它们放入 table 中，这样我就可以通过递减 "score" 对每个命中进行排序。

我还想为数据库中的每个 BLAST 命中生成一个 ORF，并将它们添加到 table 中相应的位置。

有人可以指出任何 resources/keywords 我可以搜索以了解如何操作表格数据吗？

示例：

blastn  -query genes.fasta  -subject genome.fasta  -outfmt "6 qseqid pident qcovs  evalue"

输出：

    qseqid pident qcovs evalue   
0   moaC     100.00    0.0       161.0      
1   moaC     99.38     1.0       161.0

我想从每一列中获取值并将它们用作方程式中的变量，然后在相应的行中打印该值。我将为 BLAST 使用 bash 脚本或 BioPython，因此我想将数据操作作为

的一部分

而不是这个例子的解决方案，我想看看是否有资源可以让我了解这个主题（到目前为止，我会使用电子表格程序来处理表格数据）

Answer 1

对于表格数据，我强烈推荐 pandas。

首先，您想将输出转换为 pandas DataFrame，这是一种非常适合存储表格形式数据的数据结构。

对于这个例子，我使用 tblastn 和示例文件 four_human_proteins.fasta and rhodopsin_nucs.fasta.

>>> import pandas as pd
>>> from Bio.Blast.Applications import NcbiblastnCommandline
>>> cline = NcbiblastnCommandline(cmd='/path/to/BLAST+/2.8.1/bin/tblastn', 
                          query='four_human_proteins.fasta', 
                          subject='rhodopsin_nucs.fasta', 
                          evalue='1e-10',
                          outfmt='"6 qseqid pident qcovs evalue"')
>>> print(cline)
/path/to/BLAST+/2.8.1/bin/tblastn -outfmt "6 qseqid pident qcovs evalue" -query four_human_proteins.fasta -evalue 1e-10 -subject rhodopsin_nucs.fasta

>>> blast_output = cline()[0].strip()
>>> print(blast_output) 
sp|P08100|OPSD_HUMAN    96.552  100 0.0
sp|P08100|OPSD_HUMAN    93.391  100 0.0
sp|P08100|OPSD_HUMAN    95.092  94  0.0
sp|P08100|OPSD_HUMAN    84.795  98  0.0
sp|P08100|OPSD_HUMAN    82.164  98  0.0
sp|P08100|OPSD_HUMAN    96.396  89  2.65e-67
sp|P08100|OPSD_HUMAN    92.308  89  7.50e-36
sp|P08100|OPSD_HUMAN    93.220  89  1.81e-32
sp|P08100|OPSD_HUMAN    96.296  89  6.37e-32
sp|P08100|OPSD_HUMAN    88.462  89  4.64e-12


>>> headers = ['qseqid', 'pident', 'qcovs', 'evalue']
>>> rows = [line.split() for line in blast_output.splitlines()]    
>>> df = pd.DataFrame(rows, columns=headers)
>>> print(df)
    qseqid  pident  qcovs   evalue
0   sp|P08100|OPSD_HUMAN    96.552  100     0.0
1   sp|P08100|OPSD_HUMAN    93.391  100     0.0
2   sp|P08100|OPSD_HUMAN    95.092  94  0.0
3   sp|P08100|OPSD_HUMAN    84.795  98  0.0
4   sp|P08100|OPSD_HUMAN    82.164  98  0.0
5   sp|P08100|OPSD_HUMAN    96.396  89  2.65e-67
6   sp|P08100|OPSD_HUMAN    92.308  89  7.50e-36
7   sp|P08100|OPSD_HUMAN    93.220  89  1.81e-32
8   sp|P08100|OPSD_HUMAN    96.296  89  6.37e-32
9   sp|P08100|OPSD_HUMAN    88.462  89  4.64e-12

首先我们需要告诉 pandas 哪些列包含 float。

>>> convert = {'pident': float, 
              'qcovs': float, 
              'evalue': float,
              'qseqid': str}
>>> df = df.astype(convert)

现在可以轻松地在此 DataFrame 上执行按列操作 df。定义您的评分函数并将结果添加为额外的列：

>>> df['score'] = df['qcovs'] / df['pident']  # adapt to your own needs
>>> print(df)
    qseqid  pident  qcovs   evalue  score
0   sp|P08100|OPSD_HUMAN    96.552  100.0   0.000000e+00    1.035711
1   sp|P08100|OPSD_HUMAN    93.391  100.0   0.000000e+00    1.070767
2   sp|P08100|OPSD_HUMAN    95.092  94.0    0.000000e+00    0.988516
3   sp|P08100|OPSD_HUMAN    84.795  98.0    0.000000e+00    1.155729
4   sp|P08100|OPSD_HUMAN    82.164  98.0    0.000000e+00    1.192736
5   sp|P08100|OPSD_HUMAN    96.396  89.0    2.650000e-67    0.923275
6   sp|P08100|OPSD_HUMAN    92.308  89.0    7.500000e-36    0.964163
7   sp|P08100|OPSD_HUMAN    93.220  89.0    1.810000e-32    0.954731
8   sp|P08100|OPSD_HUMAN    96.296  89.0    6.370000e-32    0.924234
9   sp|P08100|OPSD_HUMAN    88.462  89.0    4.640000e-12    1.006082

而且您可以很容易地按 score 列

对 DataFrame 进行排序

>>> df.sort_values(['score'], inplace=True)
>>> print(df)
    qseqid  pident  qcovs   evalue  score
5   sp|P08100|OPSD_HUMAN    96.396  89.0    2.650000e-67    0.923275
8   sp|P08100|OPSD_HUMAN    96.296  89.0    6.370000e-32    0.924234
7   sp|P08100|OPSD_HUMAN    93.220  89.0    1.810000e-32    0.954731
6   sp|P08100|OPSD_HUMAN    92.308  89.0    7.500000e-36    0.964163
2   sp|P08100|OPSD_HUMAN    95.092  94.0    0.000000e+00    0.988516
9   sp|P08100|OPSD_HUMAN    88.462  89.0    4.640000e-12    1.006082
0   sp|P08100|OPSD_HUMAN    96.552  100.0   0.000000e+00    1.035711
1   sp|P08100|OPSD_HUMAN    93.391  100.0   0.000000e+00    1.070767
3   sp|P08100|OPSD_HUMAN    84.795  98.0    0.000000e+00    1.155729
4   sp|P08100|OPSD_HUMAN    82.164  98.0    0.000000e+00    1.192736

学习如何使用 Bash and/or Biopython 管理表格数据（来自 BLAST+6 格式）的资源

Resources to learn how to manage tabular data (from BLAST+6 format ) using Bash and/or Biopython

python

bash

tabular

biopython

blast