從blast結果中取出每個query搜到的evalue最小的結果
阿新 • • 發佈:2018-11-21
在做多基因blast時,通常每個基因找到的匹配序列很多。這時習慣根據evalue來進行篩選,evalue較小的其相似性更高。下面提供兩種方法解決。
一 linux命令
第11列為evalue值,第一列為基因名,先根據evalue升序排列,然後根據基因名去重。預設會保留最上面的一條記錄,即evalue最小值。
二 pandas
最近在看pandas,所以拿來練手。思路也是先排序,後去重。
import pandas as pd
#將blast(-outfmt 6)輸出結果儲存到DataFrame
inp = pd.read_table( 'E:\python_test\1.blast')
inp
query | subject | identity | align_length | mismatches | gap_opens | q_start | q_end | s_start | s_end | evalue | bit_score | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | gene1 | SQ183094348 | 100 | 147 | 0 | 0 | 1 | 147 | 378 | 232 | 3 | 272 |
1 | gene1 | SQ183119192 | 100 | 66 | 0 | 0 | 1 | 66 | 82 | 147 | 2 | 122 |
2 | gene1 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 1 | 291 |
3 | gene2 | SQ183094348 | 100 | 147 | 0 | 0 | 1 | 147 | 378 | 232 | 3 | 272 |
4 | gene2 | SQ183119192 | 100 | 66 | 0 | 0 | 1 | 66 | 82 | 147 | 2 | 122 |
5 | gene2 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 1 | 291 |
6 | gene3 | SQ183094348 | 100 | 147 | 0 | 0 | 1 | 147 | 378 | 232 | 9 | 272 |
7 | gene3 | SQ183119192 | 100 | 66 | 0 | 0 | 1 | 66 | 82 | 147 | 8 | 122 |
8 | gene3 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 7 | 291 |
#取出每個query對應的evalue最低的subject
inp.sort_values(by=['query','evalue']).drop_duplicates(subset='query')
query | subject | identity | align_length | mismatches | gap_opens | q_start | q_end | s_start | s_end | evalue | bit_score | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
2 | gene1 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 1 | 291 |
5 | gene2 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 1 | 291 |
8 | gene3 | SQ182140986 | 100 | 157 | 0 | 0 | 1 | 157 | 88 | 244 | 7 | 291 |
有可能出現gene相同,evalue相同的情況,我覺得可以在加上bit_socre和align_length進行排序,這兩列為降序排列。