2015-01-10 46 views
0

我有一个包含12列的表,并希望根据第二列(sseqid)选择第一列中的项(qseqid)。这意味着第二列(sseqid)在第11列和第12列中分别以不同的值重复,分别为evaluebitscore。 ,我想获得的是那些具有最低evalue最高bitscore(当evalue可相同,列的其余部分可以被忽略,并且该数据低于下来)。Python通过使用字典比较表中的值来选择项目

所以,我做了一个简短的代码,它使用第二列作为字典的关键。我可以从列表qseqid + evalueqseqid + bitscore获得第二列的五个不同项目。

下面是代码:

#!usr/bin/python 

filename = "data.txt" 

readfile = open(filename,"r") 

d = dict() 

for i in readfile.readlines(): 
    i = i.strip() 
    i = i.split("\t") 
    d.setdefault(i[1], []).append([i[0],i[10]]) 
    d.setdefault(i[1], []).append([i[0],i[11]]) 

for x in d: 
    print(x,d[x]) 

readfile.close() 

但是,我挣扎着爬了qseqid最低的安勤和每个sseqid最高bitscore。 有没有什么好的逻辑来解决这个问题?

data.txt文件(包括标题行,并与»代表制表符)

qseqid»sseqid»pident»length»mismatch»gapopen»qstart»qend»sstart»send»evalue»bitscore 
ACLA_022040»TBB»32.71»431»258»8»39»468»24»423»2.00E-76»240 
ACLA_024600»TBB»80»435»87»0»1»435»1»435»0»729 
ACLA_031860»TBB»39.74»453»251»3»1»447»1»437»1.00E-121»357 
ACLA_046030»TBB»75.81»434»105»0»1»434»1»434»0»704 
ACLA_072490»TBB»41.7»446»245»3»4»447»3»435»2.00E-120»353 
ACLA_010400»EF1A»27.31»249»127»8»69»286»9»234»3.00E-13»61.6 
ACLA_015630»EF1A»22»491»255»17»186»602»3»439»8.00E-19»78.2 
ACLA_016510»EF1A»26.23»122»61»4»21»127»9»116»2.00E-08»46.2 
ACLA_023300»EF1A»29.31»447»249»12»48»437»3»439»2.00E-45»155 
ACLA_028450»EF1A»85.55»443»63»1»1»443»1»442»0»801 
ACLA_074730»CALM»23.13»147»101»4»6»143»2»145»7.00E-08»41.2 
ACLA_096170»CALM»29.33»150»96»4»34»179»2»145»1.00E-13»55.1 
ACLA_016630»CALM»23.9»159»106»5»58»216»4»147»5.00E-12»51.2 
ACLA_031930»RPB2»36.87»1226»633»24»121»1237»26»1219»0»734 
ACLA_065630»RPB2»65.79»1257»386»14»1»1252»4»1221»0»1691 
ACLA_082370»RPB2»27.69»1228»667»37»31»1132»35»1167»7.00E-110»365 
ACLA_061960»ACT»28.57»147»95»5»146»284»69»213»3.00E-12»57.4 
ACLA_068200»ACT»28.73»463»231»13»16»471»4»374»1.00E-53»176 
ACLA_069960»ACT»24.11»141»97»4»581»718»242»375»9.00E-09»46.2 
ACLA_095800»ACT»91.73»375»31»0»1»375»1»375»0»732 

而这里的表的内容多一点阅读的版本:

0   1   2  3  4  5  6 7  8 9  10  11 
qseqid  sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore 
ACLA_022040 TBB  32.71 431  258  8 39 468  24 423 2.00E-76  240 
ACLA_024600 TBB  80 435  87  0  1 435  1 435   0  729 
ACLA_031860 TBB  39.74 453  251  3  1 447  1 437 1.00E-121  357 
ACLA_046030 TBB  75.81 434  105  0  1 434  1 434   0  704 
ACLA_072490 TBB  41.7 446  245  3  4 447  3 435 2.00E-120  353 
ACLA_010400 EF1A 27.31 249  127  8 69 286  9 234 3.00E-13  61.6 
ACLA_015630 EF1A  22 491  255  17 186 602  3 439 8.00E-19  78.2 
ACLA_016510 EF1A 26.23 122  61  4 21 127  9 116 2.00E-08  46.2 
ACLA_023300 EF1A 29.31 447  249  12 48 437  3 439 2.00E-45  155 
ACLA_028450 EF1A 85.55 443  63  1  1 443  1 442   0  801 
ACLA_074730 CALM 23.13 147  101  4  6 143  2 145 7.00E-08  41.2 
ACLA_096170 CALM 29.33 150  96  4 34 179  2 145 1.00E-13  55.1 
ACLA_016630 CALM  23.9 159  106  5 58 216  4 147 5.00E-12  51.2 
ACLA_031930 RPB2 36.87 1226  633  24 121 1237  26 1219   0  734 
ACLA_065630 RPB2 65.79 1257  386  14  1 1252  4 1221   0  1691 
ACLA_082370 RPB2 27.69 1228  667  37 31 1132  35 1167 7.00E-110  365 
ACLA_061960 ACT  28.57 147  95  5 146 284  69 213 3.00E-12  57.4 
ACLA_068200 ACT  28.73 463  231  13 16 471  4 374 1.00E-53  176 
ACLA_069960 ACT  24.11 141  97  4 581 718 242 375 9.00E-09  46.2 
ACLA_095800 ACT  91.73 375  31  0  1 375  1 375   0  732 
+2

什么是预期的结果? – thefourtheye

回答

2
#!usr/bin/python 
import csv 

DATA = "data.txt" 

class Sequence: 
    def __init__(self, row): 
     self.qseqid =  row[0] 
     self.sseqid =  row[1] 
     self.pident = float(row[2]) 
     self.length = int(row[3]) 
     self.mismatch = int(row[4]) 
     self.gapopen = int(row[5]) 
     self.qstart = int(row[6]) 
     self.qend  = int(row[7]) 
     self.sstart = int(row[8]) 
     self.send  = int(row[9]) 
     self.evalue = float(row[10]) 
     self.bitscore = float(row[11]) 

    def __str__(self): 
     return (
      "{qseqid}\t" 
      "{sseqid}\t" 
      "{pident}\t" 
      "{length}\t" 
      "{mismatch}\t" 
      "{gapopen}\t" 
      "{qstart}\t" 
      "{qend}\t" 
      "{sstart}\t" 
      "{send}\t" 
      "{evalue}\t" 
      "{bitscore}" 
     ).format(**self.__dict__) 

def entries(fname, header_rows=1, dtype=list, **kwargs): 
    with open(fname) as inf: 
     incsv = csv.reader(inf, **kwargs) 

     # skip header rows 
     for i in range(header_rows): 
      next(incsv) 

     for row in incsv: 
      yield dtype(row) 

def main(): 
    bestseq = {} 
    for seq in entries(DATA, dtype=Sequence, delimiter="\t"): 
     # see if a sequence with the same sseqid already exists 
     prev = bestseq.get(seq.sseqid, None) 

     if (
      prev is None 
      or seq.evalue < prev.evalue 
      or (seq.evalue == prev.evalue and seq.bitscore > prev.bitscore) 
     ): 
      bestseq[seq.sseqid] = seq 

    # display selected sequences 
    keys = sorted(bestseq) 
    for key in keys: 
     print(bestseq[key]) 

if __name__ == "__main__": 
    main() 

导致在

ACLA_095800  ACT  91.73 375  31  0  1  375  1  375  0.0  732.0 
ACLA_096170  CALM 29.33 150  96  4  34  179  2  145  1e-13 55.1 
ACLA_028450  EF1A 85.55 443  63  1  1  443  1  442  0.0  801.0 
ACLA_065630  RPB2 65.79 1257 386  14  1  1252 4  1221 0.0  1691.0 
ACLA_024600  TBB  80.0 435  87  0  1  435  1  435  0.0  729.0 
0
filename = 'data.txt' 

readfile = open(filename,'r') 

d = dict() 
sseqid=[] 
lines=[] 
for i in readfile.readlines(): 
    sseqid.append(i.rsplit()[1]) 
    lines.append(i.rsplit()) 

sorted_sseqid = sorted(set(sseqid)) 

sdqDict={} 
key =None 

for sorted_ssqd in sorted_sseqid: 

    key=sorted_ssqd 
    evalue=[] 
    bitscore=[] 
    qseid=[] 
    for line in lines: 
     if key in line: 
      evalue.append(line[10]) 
      bitscore.append(line[11]) 
      qseid.append(line[0]) 
    sdqDict[key]=[qseid,evalue,bitscore] 

print sdqDict 

print 'TBB LOWEST EVALUE' + '---->' + min(sdqDict['TBB'][1]) 
##I think you can do the list manipulation below to find out the qseqid 

readfile.close() 
3

既然你是一个Python新手,我很高兴有几个例子如何手动,但为了比较,我会展示如何使用pandas库,它使表格数据简单得多。

由于您没有提供示例输出,因此我假定“每个sseqid的最低evalue和最高bitscore”指的是给定sseqid的“最低evalues中最高的bitscore”如果你想要分开,那也是微不足道的。

import pandas as pd 
df = pd.read_csv("acla1.dat", sep="\t") 
df = df.sort(["evalue", "bitscore"],ascending=[True, False]) 
df_new = df.groupby("sseqid", as_index=False).first() 

产生

>>> df_new 
    sseqid  qseqid pident length mismatch gapopen qstart qend sstart send  evalue bitscore 
0 ACT ACLA_095800 91.73  375  31  0  1 375  1 375 0.000000e+00  732.0 
1 CALM ACLA_096170 29.33  150  96  4  34 179  2 145 1.000000e-13  55.1 
2 EF1A ACLA_028450 85.55  443  63  1  1 443  1 442 0.000000e+00  801.0 
3 RPB2 ACLA_065630 65.79 1257  386  14  1 1252  4 1221 0.000000e+00 1691.0 
4 TBB ACLA_024600 80.00  435  87  0  1 435  1 435 0.000000e+00  729.0 

基本上,我们首先读取的数据文件转换成对象称为DataFrame,其是一种像Excel工作表。然后我们按evalue升序排序(以便降低evalue秒)和降低bitscore(以便更高的bitscore秒)。然后我们可以使用groupby来收集等于sseqid的数据组中的数据,并将每个组中的第一个数据作为第一个数据,这是因为排序会成为我们想要的数据。

+0

这是一个非常好的解决方案,但熊猫模块需要许多其他需要首先安装的模块。我可以在Linux上使用Python,但不幸的是,我目前正在使用Windows版本,坚持安装Six模块。无论如何,很好,谢谢。 – Karyo

2

虽然不像使用pandas库那样优雅和简洁,但很有可能在不诉诸第三方模块的情况下做你想做的事情。以下使用collections.defaultdict类来帮助创建可变长度记录列表的字典。 AttrDict类的使用是可选的,但它使访问每个基于字典的记录的字段更容易,并且不像通常所需的dict['fieldname']语法那样难看。

import csv 
from collections import defaultdict, namedtuple 
from itertools import imap 
from operator import itemgetter 

data_file_name = 'data.txt' 
DELIMITER = '\t' 
ssqeid_dict = defaultdict(list) 

# from http://stackoverflow.com/a/1144405/355230 
def multikeysort(items, columns): 
    comparers = [((itemgetter(col[1:].strip()), -1) if col.startswith('-') else 
        (itemgetter(col.strip()),  1)) for col in columns] 
    def comparer(left, right): 
     for fn, mult in comparers: 
      result = cmp(fn(left), fn(right)) 
      if result: 
       return mult * result 
     else: 
      return 0 
    return sorted(items, cmp=comparer) 

# from http://stackoverflow.com/a/15109345/355230 
class AttrDict(dict): 
    def __init__(self, *args, **kwargs): 
     super(AttrDict, self).__init__(*args, **kwargs) 
     self.__dict__ = self 

with open(data_file_name, 'rb') as data_file: 
    reader = csv.DictReader(data_file, delimiter=DELIMITER) 
    format_spec = '\t'.join([('{%s}' % field) for field in reader.fieldnames]) 
    for rec in (AttrDict(r) for r in reader): 
     # Convert the two sort fields to numeric values for proper ordering. 
     rec.evalue, rec.bitscore = map(float, (rec.evalue, rec.bitscore)) 
     ssqeid_dict[rec.sseqid].append(rec) 

for ssqeid in sorted(ssqeid_dict): 
    # Sort each group of recs with same ssqeid. The first record after sorting 
    # will be the one sought that has the lowest evalue and highest bitscore. 
    selected = multikeysort(ssqeid_dict[ssqeid], ['evalue', '-bitscore'])[0] 
    print format_spec.format(**selected) 

输出(»代表标签):

ACLA_095800» ACT» 91.73» 375» 31»  0»  1»  375» 1»  375» 0.0» 732.0 
ACLA_096170» CALM» 29.33» 150» 96»  4»  34»  179» 2»  145» 1e-13» 55.1 
ACLA_028450» EF1A» 85.55» 443» 63»  1»  1»  443» 1»  442» 0.0» 801.0 
ACLA_065630» RPB2» 65.79» 1257» 386» 14»  1»  1252» 4»  1221» 0.0» 1691.0 
ACLA_024600» TBB» 80»  435» 87»  0»  1»  435» 1»  435» 0.0» 729.0 
+0

这很好,但输出不显示qsequid。谢谢你解决我的问题。 – Karyo

+0

谢谢。几乎可笑的是,在你重新格式化你的问题之后,我是那个不了解你想要的结果的人。无论如何,我已经在这方面更新了我的答案,现在我明白什么更好。 – martineau