インデックスは、作業中の情報の短縮版です。この場合、「キー」 (@ 行の最初のコロン (':') と末尾近くの最後のスラッシュ ('/') の間のテキスト) と、何らかの値が必要になります。
この場合の「値」は 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')
このコードは、データのブロックを取得するために最初のファイルに戻ることに注意してください。ファイル間でデータが同一である場合、一致が発生したときに 2 番目のファイルからブロックをコピーできます。
また、抽出しようとしているものによっては、出力ブロックの順序を変更したり、キーが一意であることを確認したり、キーが一意ではなく、それらが一致する順序。それはあなた次第です - 私はあなたがデータで何をしているのか分かりません.