Фильтрация кадра данных последовательностей BLAST для получения в каждом кластере максимального значения pident_x - PullRequest
Новая
       13

Фильтрация кадра данных последовательностей BLAST для получения в каждом кластере максимального значения pident_x

0 голосов
/ 06 мая 2018

У меня проблема, мне нужно проанализировать следующий фрейм данных: 

    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

Пояснение к кадру данных и к выходному выводу:

  • cluster_name : - кластер, в котором присутствуют одна, две или более парных последовательностей.
  • qseqid: это имя одной последовательности
  • sseqid: это имя другой последовательности. Они делают одно сравнение qseqid против sseqid

  • pident_x: оценка после сравнения (выравнивания) этих двух последовательностей, 1 означает, что они идентичны. Когда взрыв выравнивает две последовательности, он дает мне координаты, где мои последовательности выровнены («гомологичные») в моем выравнивании, например, если у меня есть: 

             10            24

    seq1: AAATTTCCCGGGATGCGATGACGATGAAAAAATTTGG XXXXXXXXX !!!!!!!!!!!!!!! XXXXXXXXXXXXX seq2: GATGAGATCGGGATGCGATGAGGAGATAGAGATAGAG где х разница и! это совпадение, взрыв даст мне: qstart (начало первого seq): 10 qend (конец первого seq): 24 sstar (начало второй очереди): 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)

для cluster4: seq11_0035_0035 vs seq17_0042_0042 (0.97) имеет лучший результат, чем seq11_0035_0035 vs seq14_0042_0035 (0.89)

для 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

Вот мои реальные данные: данные

вот вам пример: 

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

Ответы [ 2 ]

0 голосов
/ 08 мая 2018

Вот еще один ответ, надеюсь, более эффективный.

Я сосредоточусь на главном в вашем вопросе: посмотрим, произошла ли обратная пара в паре столбцов в предыдущем ряду. В качестве примера я буду использовать 

df = pd.DataFrame({'a': [1, 1, 4, 5, 2, 5], 'b': [7, 5, 2, 1, 4, 2]})
>>> df
    a   b
0   1   7
1   1   5
2   4   2
3   5   1
4   2   4
5   5   2

Давайте найдем кортежи каждого ряда: 

this = df[['a', 'b']].apply(tuple, axis=1)
>>> this
0    (1, 7)
1    (1, 5)
2    (4, 2)
3    (5, 1)
4    (2, 4)
5    (5, 2)
dtype: object

Вот накопленная сумма появлений обратных кортежей: 

cum = pd.get_dummies(df[['b', 'a']].apply(tuple, axis=1)).cumsum()
>>> cum
(1, 5)  (2, 4)  (2, 5)  (4, 2)  (5, 1)  (7, 1)
0   0.0 0.0 0.0 0.0 1.0 1.0
1   0.0 1.0 0.0 0.0 1.0 1.0
2   1.0 1.0 0.0 0.0 1.0 1.0
3   1.0 1.0 0.0 1.0 1.0 1.0
4   1.0 1.0 1.0 1.0 1.0 1.0
5   NaN NaN NaN NaN NaN NaN

Для этого нам нужно показать, что кортеж текущей строки, если он не существует в обратном порядке, не был найден: 

this_zeros = pd.get_dummies(this)
this_zeros[:] = 0
>>> pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1)
(1.0, 5.0)  (2.0, 4.0)  (2.0, 5.0)  (4.0, 2.0)  (5.0, 1.0)  (nan, nan)  (1, 7)  (2, 2)  (5, 2)
0   0.0 0.0 0.0 0.0 1.0 0.0 0   0   0
1   0.0 1.0 0.0 0.0 1.0 0.0 0   0   0
2   1.0 1.0 0.0 0.0 1.0 0.0 0   0   0
3   1.0 1.0 0.0 1.0 1.0 0.0 0   0   0
4   1.0 1.0 1.0 1.0 1.0 0.0 0   0   0
5   1.0 1.0 1.0 1.0 1.0 1.0 0   0   0
6   NaN NaN NaN NaN NaN NaN 0   0   0

Теперь мы можем использовать этот DataFrame для поиска текущего кортежа: 

keep = pd.concat([cum, this_zeros[this_zeros.columns.difference(cum.columns)]], axis=1).lookup(df.index, this)

Мы должны сохранить исходный фрейм данных 

>>> df[keep.astype(bool)]
    a   b
0   1   7
1   1   5
2   4   2
5   5   2
0 голосов
/ 06 мая 2018

Если вы создаете кортеж из столбцов, а затем выполняете кумулятивную сумму, вы можете проверить, появляется ли обратная пара в кумулятивной сумме: 

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) создает обратные кортежи и кумулятивно суммирует их в кортежах

  • Остальное проверяет, входит ли одно в другое.
...