python - 过滤 BLAST 序列的数据帧以在每个集群中获得最大 pident_x

Question

我有一个问题，我需要解析以下数据框:

    cluster_name    qseqid  sseqid  pident_x    qstart  qend    sstar   send
2   1   seq1_0035_0035  seq13_0042_0035 0.73    42  133 46  189
3   1   seq1_0035_0035  seq13_0042_0035 0.73    146 283 287 389
4   1   seq1_0035_0035  seq13_0042_0035 0.73    301 478 402 503
5   1   seq13_0042_0035 seq1_0035_0035  0.73    46  189 42  133
6   1   seq13_0042_0035 seq1_0035_0035  0.73    287 389 146 283
7   1   seq13_0042_0035 seq1_0035_0035  0.73    402 503 301 478
8   2   seq4_0042_0035  seq2_0035_0035  0.71    256 789 125 678
9   2   seq4_0042_0035  seq2_0035_0035  0.71    802 1056    706 985
10  2   seq4_0042_0035  seq7_0035_0042  0.83    123 745 156 723
12  4   seq11_0035_0035 seq14_0042_0035 0.89    145 647 236 921
13  4   seq11_0035_0035 seq17_0042_0042 0.97    148 623 241 1002
14  5   seq17_0035_0042 seq17_0042_0042 0.94    188 643 179 746

dataframe 和 blast 输出的解释:

• cluster_name :是存在一个、两个或更多配对序列的簇。
• qseqid : 是一个序列的名称
• sseqid : 是另一个序列的名称。这些做一个比较 qseqid 与 sseqid

• pident_x : 是这两个序列比对(比对)后的得分，1表示相同。当 blast 比对两个序列时，它会给出我的序列在我的比对中对齐(“同源”)的坐标，例如，如果我有:

           10            24
  

  seq1 : AAATTTCCCGGGATGCGATGACGATGAAAAAATTTGG xxxxxxxxx!!!!!!!!!!!!!!!!xxxxxxxxxxxxxseq2:GATGAGATCGGGATGCGATGAGGAGATAGAGATAGAG其中 x 和 !是一场比赛，爆炸会给我:qstart(第一个序列的开始):10qend(第一个序列结束):24sstar(第二个序列的开始):10发送(第二个序列结束):24

注意:这是一个示例，但不一定从 0 开始。

我真正想要的只是在每个集群中获得最大的 pident_x 但问题是如你所见我可以有相反的序列(如果你看一下 2,3,4 和 5,6,7它们是相同的但相反)，我需要做的是只保留一个，例如第 2,3 和 4 行，因为 blast 将比较每个序列，甚至是倒数序列。

那么输出将是:

cluster_name    qseqid  sseqid  pident_x    qstart  qend    sstar   send
    2   1   seq1_0035_0035  seq13_0042_0035 0.73    42  133 46  189
    3   1   seq1_0035_0035  seq13_0042_0035 0.73    146 283 287 389
    4   1   seq1_0035_0035  seq13_0042_0035 0.73    301 478 402 503
    10  2   seq4_0042_0035  seq7_0035_0042  0.83    123 745 156 723
    13  4   seq11_0035_0035 seq17_0042_0042 0.97    148 623 241 1002
    14  5   seq17_0035_0042 seq17_0042_0042 0.94    188 643 179 746

确实:对于 cluster1: seq1_0035_0035 vs seq13_0042_0035他的反转了seq13_0042_0035 seq1_0035_0035但我只保留第一个。

对于 cluster2: seq4_0042_0035 vs seq7_0035_0042 (0.83)具有比 seq4_0042_0035 vs seq2_0035_0035 (0.71) 更好的 pident 分数

对于 cluster4: seq11_0035_0035 vs seq17_0042_0042 (0.97)具有比 seq11_0035_0035 vs seq14_0042_0035 (0.89) 更好的 pident 分数

对于 custer5:只有一个配对序列 seq17_0035_0042 vs seq17_0042_0042 (0.94) ，那我留着这个

我真的不知道如何做这样的事情，有人有想法吗？

部分补充:

这是我从这个小数据集使用的脚本(与我上面的示例相同):smalldata

blast=pd.read_csv("match.csv",header=0)

#blast=blast.drop(columns=[ "length", "mismatch", "gapopen", "evalue", "bitscore","pident"])

#Cluster Dataframe
cluster=pd.read_csv("cluster_test.csv",header=0)
cluster.columns = ["cluster_name", "seq_names"]

#Distance mean dataframe
dist=pd.read_csv("fnode.csv",header=0)
dist.columns = ["qseqid", "sseqid","pident","coverage"]
dist=dist.drop(columns=["coverage"])

#Including cluster information and distance mean information into one dataframe:
data = cluster.merge(dist, left_on='seq_names', right_on='qseqid')

#Adding for each two remaining dataframe a concatened colomn
data["name_concatened"] = data["qseqid"].map(str) + data["sseqid"]
blast["name_concatened"] = blast["qseqid"].map(str) + blast["sseqid"]
#We do not need these columns anymore
blast=blast.drop(columns=[ "qseqid","sseqid"])

#Including cluster information + distance mean information  + coordinate sequences from blast into one dataframe:
data = data.merge(blast, left_on='name_concatened', right_on='name_concatened')
data=data.drop(columns=[ "seq_names","name_concatened","pident_y"])

print(data)

this = data[["qseqid", "sseqid"]].apply(tuple, axis=1)
cum = pd.get_dummies(data[["sseqid", 'qseqid']].apply(tuple, axis=1)).cumsum()

this_zeros = pd.get_dummies(this)
this_zeros[:] = 0
pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1)
keep = pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1).lookup(data.index, this)

data=data[keep.astype(bool)]

print(data)

但是正如你在这里看到的，我只得到:

  cluster_name           qseqid          sseqid  pident_x  qstart  qend  \
4             1  seq13_0042_0035  seq1_0035_0035      0.73      46   189   
5             1  seq13_0042_0035  seq1_0035_0035      0.73     287   389   
6             1  seq13_0042_0035  seq1_0035_0035      0.73     402   503   

   sstar  send  
4     42   133  
5    146   283  
6    301   478   

我应该得到:

cluster_name    qseqid  sseqid  pident_x    qstart  qend    sstar   send
        2   1   seq1_0035_0035  seq13_0042_0035 0.73    42  133 46  189
        3   1   seq1_0035_0035  seq13_0042_0035 0.73    146 283 287 389
        4   1   seq1_0035_0035  seq13_0042_0035 0.73    301 478 402 503
        10  2   seq4_0042_0035  seq7_0035_0042  0.83    123 745 156 723
        13  4   seq11_0035_0035 seq17_0042_0042 0.97    148 623 241 1002
        14  5   seq17_0035_0042 seq17_0042_0042 0.94    188 643 179 746

这是我的真实数据:datas

举个例子:

df = pd.DataFrame({'a': [1, 1, 4, 5, 2, 5], 'b': [7, 5, 2, 1, 4, 2]})
this = df[['a', 'b']].apply(tuple, axis=1)
cum = pd.get_dummies(df[['b', 'a']].apply(tuple, axis=1)).cumsum()
this_zeros = pd.get_dummies(this)
this_zeros[:] = 0
pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1)
keep = pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1).lookup(df.index, this)
df=df[keep.astype(bool)]
print(df)

我的结果:

 a  b
3  5  1
4  2  4

Answer 1

如果您从列中创建一个元组，然后执行累加和，您可以检查反向对是否已经出现在累加和中:

df[~pd.DataFrame({
    'tup': df[['sseqid', 'qseqid']].apply(tuple, axis=1), 
    'inv_tups': df[['qseqid', 'sseqid']].apply(lambda t: (tuple(t), ), axis=1).cumsum().shift(1)}
).apply(lambda r: isinstance(r.inv_tups, tuple) and r.tup in r.inv_tups, axis=1)]

• df[['sseqid', 'qseqid']].apply(tuple, axis=1) 从列中创建元组。

• df[['qseqid', 'sseqid']].apply(lambda t: (tuple(t), ), axis=1).cumsum().shift(1) 创建逆元组，并将它们累加到元组中

• 其余的检查一个是否包含在另一个中。