分染色体并行计算中科学分配任务

在生物信息分析过程中,为了加快运行速度,我们经常会采用分染色体计算的方式,即scatter/gather的思想。比如SNP calling, SNP phasing及imputation等。但怎样划分染色体,划分几个部分比较合适呢?

一般情况下,程序运行的时间和染色体长度成正比,而组装出来的染色体一般长度差异很大。在计算资源充足的情况下,一条染色体一个任务通常是最快的。如果任务存在排队,包含长染色体的任务如果排在后面,那需要的时间反而更多了。之前供职的公司,流程采用随机分成染色体数目相等的N份,但存在一些问题:1)N为多少合适?多了如果任务有排队,一条染色体一个任务不一定快,还比较浪费计算资源。分太少又需要更多计算时间。2)随机分配,很有可能将两条较长的染色体分在一个任务中。

很多人想到了,既然程序运行时间和染色体长度成正比,那保证每份任务所分配的染色体总长相等不就可以了吗?确实是这样,这里有一个前提,即染色体不能打断,因为打断后,有些计算可能出问题,比如断点处存在变异,比如计算ROH、IBD等。

这里提供一种比较科学的解决办法:将该问题视为简单的一维装箱问题,箱子的大小即最长染色体的长度,利用贪婪算法解决这个装箱问题。如果最长的染色体计算时间是最长的,这样其它任务都不会超过这个时间,而且保证了分配的任务数量相对少(贪婪算法的结果很有可能不是全局最优解,但是好在算法简单,相对最优的结果也可以接受)。python实现如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def SplitFai(fai, outdir, size=0):
chrs = []
flist = []
for line in open(fai):
info = line.strip().split()[0:2]
chrs += [[info[0], int(info[1])]]
chrs = sorted(chrs, key=lambda a: a[1], reverse=True)
if not size:
size = chrs[0][1]
elif size < chrs[0][1]:
sys.stderr.write('Error: size is short than the longest chromosome!')
sys.exit(0)
count = 0
while chrs:
fo = open(outdir+'/gemone.split_'+str(count)+'.list', 'w')
flist.append(outdir+'/gemone.split_'+str(count)+'.list')
lenC = 0
chrs_new = []
for each in chrs:
if size - lenC >= each[1]:
fo.write(each[0] + '\n')
lenC += each[1]
else:
chrs_new += [each]
chrs = chrs_new
count += 1
return flist

该函数有3个参数,第1个是参考序列的index文件,用samtools faidx ref.fa生成,该文件第一列为染色体名,第二列为染色体长度。第2个参数为输出目录,第3个参数为每份染色体总长度的上限,即箱子大小,默认为最长染色体长度。 注意,如果自己指定箱子大小,箱子大小必须超过最长染色体长度。 输出文件为每份包含的染色体ID。

该函数的思路:将所有染色体按照长度从大到小排序。装箱的时候先尽量装大的,直到装不下再换新箱子。

如果我的博客能帮你节约时间,我会非常高兴你请我喝杯咖啡。