索引只是您正在使用的信息的縮短版本。在這種情況下，您需要「鍵」 - @ -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