双循环需要时间

Question

我有一个脚本需要花费很多时间,并且在2天之后无法完成...我将1个文件解析为2个字典,如下所示:

gfftree = {'chr1':[(gene_id, gstart, gend),...], 'chr2':[(gene_id, gstart, gend),...],...}
TElocation = {'chr1':[(TE_id, TEstart, TEend),...], 'chr2':[(TE_id, TEstart, TEend),...],...}

.

- 目的是找到TE_id,其TEstart或TEend或两者都位于每个chr(密钥)中的gene_id'gstart和gend之间.

以上内容应更改为"查找TE_id,其范围(TEstart,TEend)与任何gene_id的范围重叠(gstart,gend)"

这是我的代码:

TE_in_TSS = []
for TErange in TElocation[chromosome]:
    TE_id, TEstart, TEend = TErange
    for item in gfftree[chromosome]:
        gene, gstart, gend = item       
        if len(list(set(range(int(gstart),int(gend)+1)) & set(range(int(TEstart),int(TEend)+1)))) > 0:
            TE_in_TSS.append((gene, TE_id, TEstart, TEend))
        else:
            pass

到目前为止,我确信这个循环对于小数据很好,但是当涉及到800,000 TE_id和4,000 gene_id等更大的循环时,需要时间......而且我不知道它是否可以完成......

Answer 1

OP 方法是O(n*m)，其中n是基因数量，m是 TE 数量。这种方法不是像OP中那样针对每个TE测试每个基因，而是利用基因和TE的有序性质以及指定的匹配规则，仅查看每个基因和TE一次，除了所描述的基因前瞻之外在3.下面。这种方法的O(n + m)前提是平均基因前瞻相对于 来说较小n。每个基因和 TE 被访问的序列描述如下：

1. 当我们完成针对当前基因的当前 TE 测试后，我们得到下一个 TE。
2. 当当前 TE 的起始位置超过当前基因的结束位置时，我们得到下一个基因，直到不是。
3. 如果我们找到匹配的 TE/基因对，我们将针对当前 TE 测试每个连续的基因，直到没有匹配，保持当前基因不变。

def get_TE_in_TSS(genes, TEs):
    TE_in_TSS = []
    gene_pos, TE_pos = 0, 0
    gene_count, TE_count = len(genes), len(TEs)
    while gene_pos < gene_count:
        while (TE_pos < TE_count) and (TEs[TE_pos][1] <= genes[gene_pos][2]):
            match_gene_pos = gene_pos
            while (match_gene_pos < gene_count) and (TEs[TE_pos][2] >= genes[match_gene_pos][1]):
                TE_in_TSS.append((genes[match_gene_pos][0], TEs[TE_pos][0],
                                  TEs[TE_pos][1], TEs[TE_pos][2]))
                match_gene_pos += 1 # look ahead to see if this TE matches the next gene
            TE_pos += 1
        gene_pos += 1
    return TE_in_TSS

OP 报告的性能：

1 second (compared to 2 days + for OP code) for 801,948 TEs, 6,007 genes

测试数据：

genes = (('HTR3A', 7, 9), ('ADAMTSL4', 10,100), ('THSD4',2000, 2800), ('PAPLN', 2850, 3000))
TEs = (('a', 10, 11), ('b', 13, 17), ('c', 50, 2500), ('d', 2550, 2700),
       ('e', 2800, 2900), ('f', 9999, 9999)) 
TE_in_TSS = get_TE_in_TSS(genes, TEs)
print(TE_in_TSS)

输出：

[('ADAMTSL4', 'a', 10, 11), ('ADAMTSL4', 'b', 13, 17), ('ADAMTSL4', 'c', 50, 2500), 
 ('THSD4', 'c', 50, 2500), ('THSD4', 'd', 2550, 2700), ('THSD4', 'e', 2800, 2900), 
 ('PAPLN', 'e', 2800, 2900)]

请注意，本文的前 9 条评论提到了一种更有效的O(n * m)方法，该方法因明确的规范而变得过时。