索引只是您正在使用的信息的缩短版本。在这种情况下,您需要“键” - @ -line的第一个冒号(':')和最后的最后一个斜杠('/')之间的文本以及某种类型的值。由于在这种情况下的“值”是4行块的全部内容,并且由于我们的索引将为每个块存储单独的条目,所以如果我们使用了,我们将把整个文件存储在内存中索引中的实际值。
取而代之,让我们使用4行块开头的文件位置。这样,您可以移动到该文件位置,打印4行,然后停止。总成本是存储整数文件位置所需的4或8或多个字节,而不是实际基因组数据的许多字节。
下面是一些代码,完成这项工作,但也做了大量的验证和检查。你可能想扔掉你不用的东西。
import sys
def build_index(path):
index = {}
for key, pos, data in parse_fastq(path):
if key not in index:
# Don't overwrite duplicates- use first occurrence.
index[key] = pos
return index
def error(s):
sys.stderr.write(s + "\n")
def extract_key(s):
# This much is fairly constant:
assert(s.startswith('@'))
(machine_name, rest) = s.split(':', 1)
# Per wikipedia, this changes in different variants of FASTQ format:
(key, rest) = rest.split('/', 1)
return key
def parse_fastq(path):
"""
Parse the 4-line FASTQ groups in path.
Validate the contents, somewhat.
"""
f = open(path)
i = 0
# Note: iterating a file is incompatible with fh.tell(). Fake it.
pos = offset = 0
for line in f:
offset += len(line)
lx = i % 4
i += 1
if lx == 0: # @machine: key
key = extract_key(line)
len1 = len2 = 0
data = [ line ]
elif lx == 1:
data.append(line)
len1 = len(line)
elif lx == 2: # +machine: key or something
assert(line.startswith('+'))
data.append(line)
else: # lx == 3 : quality data
data.append(line)
len2 = len(line)
if len2 != len1:
error("Data length mismatch at line "
+ str(i-2)
+ " (len: " + str(len1) + ") and line "
+ str(i)
+ " (len: " + str(len2) + ")\n")
#print "Yielding @%i: %s" % (pos, key)
yield key, pos, data
pos = offset
if i % 4 != 0:
error("EOF encountered in mid-record at line " + str(i));
def match_records(path, index):
results = []
for key, pos, d in parse_fastq(path):
if key in index:
# found a match!
results.append(key)
return results
def write_matches(inpath, matches, outpath):
rf = open(inpath)
wf = open(outpath, 'w')
for m in matches:
rf.seek(m)
wf.write(rf.readline())
wf.write(rf.readline())
wf.write(rf.readline())
wf.write(rf.readline())
rf.close()
wf.close()
#import pdb; pdb.set_trace()
index = build_index('afile.fastq')
matches = match_records('bfile.fastq', index)
posns = [ index[k] for k in matches ]
write_matches('afile.fastq', posns, 'outfile.fastq')
请注意,此代码将返回到第一个文件以获取数据块。如果您的数据在文件之间是相同的,那么您可以在发生匹配时从第二个文件复制该块。
还要注意,根据您要提取的内容,您可能需要更改输出块的顺序,并且您可能需要确保密钥是唯一的,或者可能确保密钥不是唯一的但是按照它们匹配的顺序重复。这取决于你 - 我不确定你在处理数据。
>非常感谢。对于像我这样的初学者来说,这是一件好事,但查看代码对我的学习很有帮助。我将afile和bfile添加为'sys.argv [1]'和'sys.argv [2]'作为命令行运行(通常从未遇到过这样的问题)。但是,当我尝试运行它时,它会显示'> /some/path/code.py(92)()'' - > index = build_index(afile)''(Pdb)'并且仍然停留在那里。它不会给出错误,它只是保持卡住... –
biohazard
对不起!我在代码中调用了调试器。带有“import pdb; pdb.set_trace()”的行导致代码停止并调用调试器。注释掉或删除该行,它应该贯穿始终。 (或者,使用“n”表示下一步,“s”表示步入,“p expr”表示打印表达式,以便在运行时观察代码。) –
现在,我收到以下警告并且没有输出。 。我很抱歉打扰你: 'Traceback(最近调用最后一次):' '文件“py_fetch_pair.py”,第94行,在' 'index = build_index(afile)' 'File (路径)中的key,pos,data中的“py_fetch_pair.py”,第12行,build_index'':' '文件“py_fetch_pair.py”,第32行,parse_fastq'中 'f = open(path) ' 'TypeError:强制为Unicode:需要字符串或缓冲区,找到文件' –
biohazard