多文件比较

Question

我想比较多个文件(15-20)，它们是gzipped的，并从它们中还原出常见的行。但这并不是那么简单。在某些列中是精确的行，而且我也希望它们在它们存在的文件中有多少个信息。如果1，这一行是唯一的文件，等等，也会很好地保存这些文件的名称。

每个文件看起来如下所示：

##SAMPLE=<ID=NormalID,Description="Cancer-paired normal sample. Sample ID 'NORMAL'">
##SAMPLE=<ID=CancerID,Description="Cancer sample. Sample ID 'TUMOR'">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NormalID_NORMAL CancerID_TUMOR
chrX    136109567       .       C       CT      .       PASS    IC=8;IHP=8;NT=ref;QSI=35;QSI_NT=35;RC=7;RU=T;SGT=ref->het;SOMATIC;TQSI=1;TQSI_NT=1;phastCons;CSQ=T|ENSG00000165370|ENST00000298110|Transcript|5KB_downstream_variant|||||||||YES|GPR101|||||        DP:DP2:TAR:TIR:TOR:DP50:FDP50:SUBDP50   23:23:21,21:0,0:2,2:21.59:0.33:0.00   33:33:16,16:13,13:4,4:33.38:0.90:0.00
chrX    150462334       .       T       TA      .       PASS    IC=2;IHP=2;NT=ref;QSI=56;QSI_NT=56;RC=1;RU=A;SGT=ref->het;SOMATIC;TQSI=2;TQSI_NT=2;CSQ=A||||intergenic_variant||||||||||||||| DP:DP2:TAR:TIR:TOR:DP50:FDP50:SUBDP50 30:30:30,30:0,0:0,0:31.99:0.00:0.00     37:37:15,17:16,16:6,5:36.7:0.31:0.00

文件是以选项卡分隔的。如果行以#开头，则忽略这一行。我们只对那些不感兴趣的问题感兴趣。以0为基础的python坐标，我们对0,1,2,3,4字段感兴趣。它们必须在被报告为普通文件的文件之间匹配。但是，我们仍然需要保存有关库仑/字段其余部分的信息，以便可以将它们写入输出文件。

现在我有以下代码：

import gzip
filenames = ['a','b','c']
files = [gzip.open(name) for name in filenames]

sets = [set(line.strip() for line in file if not line.startswith('#')) for file in files]
common = set.intersection(*sets)
for file in files: file.close()
print common

在我的当前代码中，我不知道如何正确实现if not line.startswith() (哪个位置？)，以及如何指定应该匹配的列。更不用说，我不知道如何获得行，例如，在6个文件中，或在总共15个文件中的10个中。对此有什么帮助吗？

Answer 1

收集字典中的行，其中的字段使它们与键相似：

from collections import defaultdict
d = defaultdict(list)

def process(filename, line):
    if line[0] == '#':
        return

    fields = line.split('\t')
    key = tuple(fields[0:5]) # Fields that makes lines similar/same
    d[key].append((filename, line))

for filename in filenames:
    with gzip.open(filename) as fh:
        for line in fh:
            process(filename, line.strip())

现在，您有了一个包含文件名行元组列表的字典。现在，您可以打印所有出现超过10次的线条：

for l in d.values():
   if len(l) < 10: continue

   print 'Same key found %d times:' % len(l)

   for filename, line in l:
       print '%s: %s' % (filename, line)