#!/usr/bin/env python
# -*- coding: utf-8 -*-
infile2 = open('genemark.gff3', 'r')
infile1 = set(line1.strip() for line1 in open('1.txt', 'r'))
for line in infile2:
line = line.strip().split()
if line[2] == 'gene':
chr, start, end = line[0], int(line[3]), int(line[4])
for line1 in infile1:
line1 = line1.split()
chr1, start1, end1 = line1[1], int(line1[2]), int(line1[3])
if chr1 == chr:
if start1 < start < end1:
print line1[0], line[-1]
if start1 < end < end1:
print line1[0], line[-1]
if start1 > start and end > end1:
print line1[0], line[-1]
genemark.gff3
格式类似下边:
chr1D GeneMark.hmm gene 2705930 2711118 . + . ID=1903228_g;Name=1903228_g
chr1D GeneMark.hmm mRNA 2705930 2711118 . + . ID=1903228_t;Name=1903228_t;Parent=1903228_g
1.txt
:
UN011157 chr1D 2705329 2706342 98.4 95.0 972 30 21 0
UN003843 chr1D 2705681 2721144 61.4 97.4 633 12 5 0
附上原始文件的百度云链接,希望感兴趣的参考
点击下载 密码 enu8
综合楼下各位朋友的答案,现推荐两种
第一种 根据 @ferstar @用筹兮用严 的答案,即并行版
#!/usr/bin/env python
# encoding: utf-8
from collections import defaultdict
from multiprocessing import Pool, cpu_count
from functools import partial
def find_sth(f2, f1=None):
start, end = int(f2[3]), int(f2[4])
for uno1, start1, end1 in f1[f2[0]]:
if (start1 <= start and start <= end1) or (start1 <= end and end <= end1) or (start1 >= start and end >= end1):
with open("out.txt", "a") as fh:
fh.write(uno1 + "\t" + f2[-1] + "\n")
#print(uno1, f2[-1])
def main():
with open('1.txt', 'r') as f1:
infile1 = defaultdict(set)
for uno1, chr1, start1, end1, *others in map(str.split, f1):
infile1[chr1].add((uno1, int(start1), int(end1)))
with open('genemark.gff3', 'r') as f2:
infile2 = [x for x in map(str.split, f2) if x[2] == 'gene']
pool = Pool(cpu_count())
pool.map(partial(find_sth, f1=infile1), infile2)
pool.close()
pool.join()
if __name__ == "__main__":
main()
第二种 @citaret 他的版本(单核版),对单核来说,不逊于上述代码。但是两者结果稍有不同,并行版结果更全(这里少了73条,出在判断条件的边界问题,由于对intervaltree熟悉,怎么改还不知道),现在边界问题已修改,两种代码结果完全一样,perfect!
如下
from collections import defaultdict
from intervaltree import Interval, IntervalTree
with open('1.txt') as f:
d1 = defaultdict(list)
xs = map(lambda x: x.strip().split(), f)
for x in xs:
y = (x[0], int(x[2]), int(x[3]))
d1[x[1]].append(y)
for k, v in d1.items():
d1[k] = IntervalTree(Interval(s, e, u) for u, s, e in v)
with open('genemark.gff3') as f:
for line in f:
line = line.strip().split()
if line[2] == 'gene':
chr, start, end = line[0], int(line[3]), int(line[4])
for start1, end1, un1 in d1[chr][start-1:end+1]:
print(un1, line[-1])
阿神2017-04-18 09:18:12
ここに 2 つの提案があります:
コードのネストが深すぎます。関数内でできるだけ早く返すことにより、ネスト レベルを下げることができます。同様に、ループ内で continue を使用してネスト レベルを下げることができます。
パフォーマンスについて
ループを通過するたびに file1 の行を分割するのは非常に賢明ではありません
以下は私が変更したコードです
リーリーPHPz2017-04-18 09:18:12
時間とスペースを交換して、genemark.gff3 のリストと 1.txt の辞書をそれぞれ構築します。具体的な実装:
リーリー修正バージョン v2 では、内部ループの int() が削除され、出力が簡素化されています。
リーリーv3: 質問の意味を注意深く検討した結果、メインループはフラグメントと交差するセット内のすべてのフラグメントを見つけることであることがわかりました。まずこのセットを見てみましょう。
リーリー各コレクション内のフラグメントの数は 6000 ~ 10000 であり、トラバーサルは非効率的です。そのため、フラグメントと交差するすべてのフラグメントを迅速に取得するには、intervaltree の使用を検討してください。
リーリー時間テストの結果: インターバルツリーの構築には 10 秒かかりましたが、交差処理の速度は約 100 倍向上しました。
間隔リファレンス https://pypi.python.org/pypi/...
伊谢尔伦2017-04-18 09:18:12
非常に興味深いものを見つけました。皆さんは非常に好意的な反応を示しましたが、実際の結果については、退屈でちょっとテストしてみただけです。
に二重化しました。 リーリー質問で提供されたサンプルテキストは 2 行しかないため、
1.txt
とgenemark.gff3
をそれぞれ4000
行
回答した階数で並べ替えます。たとえば、質問者のコードは hi.py
、次に 1 階の回答者のコードは hi1.py
というようになります。
概要
私の計画 - 4 階のコードに小さな変更を並行して加えます
私が書いたものに問題がありました。@yongchixiyongyan が正しい並列コードを更新しました。私のコードは変更されません。そのため、後で見るクラスメートが参照できるようになります直接採点 (Python3)
リーリー
次に、運用効率を見てみましょうリーリー
時間の点ではかなり遅いようです (4000 行のデータはわずか数百 KB です)。質問者は、処理されるデータが大きくなるほど、より明らかになります。並列処理による効率の利点追伸: 質問者が処理するデータの実際のサイズは、MB または GB レベルに達するはずだと推測しています。このレベルでの並列処理が適切な方法です。
ソースデータと結果のアドレスリンク: http://pan.baidu.com/s/1hrSZQuS パスワード: u93n