最近师兄给了小任务,算一个P值。任务详情是这样的:第一步,有基因集A,23个元素,基因集B,451个元素,共有23*451=9922个组合(当然你要考虑去重),看在PPI数据库(145万多行数据)中出现的组合总个数(观察值);第二步,然后再以hg19中基因为背景基因集,放回随机抽取23个基因构成基因集C(记得要先把451个基因在hg19基因集中去掉),与集合B组合,看在PPI数据库中出现的组合总个数(期望值)。第二步重复一千万次,计算期望值大于观察值的总次数。
师兄给了一个perl脚本,说有bug,让我调试。当我看到那个脚本使用就硬生生循环一千万次时候,我决定自己造轮子。
多进程肯定是必须的。
在第一个版本的,多进程是这样的。
</pre><pre name="code" class="python">def run(EXP_BG,PPI,BG_gene2,nu):
EXP=0
tmp_gene=[]
#EXP_BG 是hg19中基因
gene1=random.sample(EXP_BG,23)
gene2=BG_gene2#背景不变=gene[in2]
#gene2=random.sample(EXP_BG,n2)
for i in gene1:
for j in gene2:
if i == j :
<span style="white-space:pre"> </span>continue
if [i,j] in tmp_gene or [j,i] in tmp_gene:
continue
tmp_gene.append([i,j])
for i in tmp_gene:
##cycle 10373
for j in PPI:
<span style="white-space:pre"> </span>if i in PPI or [i[1],i[0]] in PPI:
<span style="white-space:pre"> </span>EXP+=1
output_dict[nu]=EXP
time1=int(time.time())#start
manager=Manager()
output_dict = manager.dict()
pool=Pool(processes=cpu)
for nu in range(10000000):
#gene1=random.sample(<span style="font-family: Arial, Helvetica, sans-serif;">EXP_BG,PPI,BG_gene2,23)</span>
pool.apply_async(run,(nu,))
pool.close()
pool.join()
</pre>当然这种方式也是能充分利用CPU的,当时没有考虑效率,算了12小时,依然没有出结果,师兄说正常半天就完事了。接着就改进吧。</p><p><pre name="code" class="python">def run(EXP_BG,PPI,BG_gene2,nu):
EXP=0
tmp_gene=[]
#EXP_BG 是hg19中基因
gene1=random.sample(EXP_BG,23)
gene2=BG_gene2#背景不变=gene[in2]
#gene2=random.sample(EXP_BG,n2)
tmp_gene=[[i,j] for i in gene1 for j in gene2 if [i,j] in PPI or [j,i] in PPI] #cycle 10373
output_dict[nu]=len(tmp_gene)
time1=int(time.time())#start
manager=Manager()
output_dict = manager.dict()
pool=Pool(processes=cpu)
for nu in range(10000000):
#gene1=random.sample(EXP_BG,23)
pool.apply_async(run,(<span style="font-family: Arial, Helvetica, sans-serif;">EXP_BG,PPI,BG_gene2,23</span><span style="font-family: Arial, Helvetica, sans-serif;">)</span>
pool.close()
pool.join()
看起来简省很多,但是实际工作量也并没有减少,一千万次*9922次循环,好吧,还是给跪了,循环还不是问题,最要的是那个in判断语句,很是拖慢速度。其实整个run函数做的就是两个集合取交集,然后大神告诉我,A&B 比[ i for i in A if i in B]快的飞起。然后我就用上了Python的set函数
def run(EXP_BG,PPI,gene2,nu):
EXP=0
tmp_gene=[]
gene1=random.sample(EXP_BG,23)
gene2=BG_gene2#背景不变=gene[in2]
tmp_gene=["__".join(sorted([i,j])) for i in gene1 for j in gene2 ]
tmp_gene=set(tmp_gene)
output_dict[nu]=len(tmp_gene&PPI)
##交集
然后就飞起,效率按理说应该很快啊,也确实是这样。但是在CPU使用上,却更坑爹了。我在想是不是每次向每个进程传递的参数和获取参数这个耽误时间了,毕竟,有两个参数是集合,一个元素个数是2万+,一个是145万+。当我尝试用manager.list()来给进程提供这些集合的时候,程序跑得飞快,3个小时左右跑完,验证的我的猜想。
<pre name="code" class="python"><pre name="code" class="python">def run(nu):
EXP=0
tmp_gene=[]
gene1=random.sample(EXP_BG,23)
gene2=BG_gene2#背景不变=gene[in2]
tmp_gene=["__".join(sorted([i,j])) for i in gene1 for j in gene2 ]
tmp_gene=set(tmp_gene)
output_dict[nu]=len(tmp_gene&PPI)
time1=int(time.time())#start
manager=Manager()
<pre name="code" class="python">PPI_list=manager.list()
PPI_list=PPI
BG_gene2_list=manager.list()
BG_gene2_list=BG_gene2
EXP_BG_list=manager.list()
EXP_BG_list=EXP_BG
output_dict = manager.dict() pool=Pool(processes=cpu)for nu in range(10000000):#gene1=random.sample(
EXP_BG,PPI,BG_gene2,23)pool.apply_async(run,(nu,))pool.close()pool.join()