回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>编辑以澄清输入/输出。我认为不管怎样,这都会有点慢,但是到目前为止,我还没有真正考虑过python脚本的速度问题,我正试图找出加快这些操作的方法。在</p>
<p>我输入的是基因组序列的词典。我目前正在研究两个基因组,萌芽酵母基因组(磁盘上11.5 MB)和人类基因组(磁盘上2.8 GB)。这些词典的形式如下:</p>
<pre><code>seq_d = { 'chr1' : 'ATCGCTCGCTGCTCGCT', 'chr2' : 'CGATCAGTCATGCATGCAT',
'chr3' : 'ACTCATCATCATCATACTGGC' }
</code></pre>
<p>我想在基因组的两条链中找到核苷酸的所有单碱基实例。其中<code>+</code>链指上述字典中的序列,<code>-</code>链是序列的反向补码。我的输出是一个嵌套字典,其中顶层<code>keys</code>是<code>+</code>或{<cd2>},嵌套的{<cd3>}是染色体名称,值是0索引位置的列表:</p>
^{pr2}$
<p><code>test_d</code>定义了一组要在脚本后面的大型Illumina排序数据集中检查的位置集。在</p>
<p>我的第一次尝试使用<code>enumerate</code>,和迭代。在</p>
<pre><code>import time
import numpy as np
rev_comps = { 'A' : 'T', 'T' : 'A', 'G' : 'C', 'C' : 'G', 'N' : 'N'}
test_d = { '+' : {}, '-' : {}}
nts = 'T'
s = time.time()
for chrom in seq_d:
plus_pos, minus_pos = [], []
chrom_seq = seq_d[chrom]
for pos, nt in enumerate(chrom_seq):
if nt in nts:
plus_pos.append(pos)
if rev_comps[nt] in nts:
minus_pos.append(pos)
test_d['+'][chrom] = plus_pos
test_d['-'][chrom] = minus_pos
e = time.time()
print 'The serial version took {} minutes...'.format((e-s)/60)
</code></pre>
<p>酵母产量:</p>
<pre><code>The serial version took 0.0455190300941 minutes...
</code></pre>
<p>人类的输出:</p>
<pre><code>The serial version took 10.1694028815 minutes...
</code></pre>
<p>我尝试使用numpy数组而不是<code>enumerate()</code>和迭代:</p>
<pre><code>s = time.time()
for chrom in seq_d:
chrom_seq = np.array(list(seq_d[chrom]))
nts = list(nts)
rev_nts = [rev_comps[nt] for nt in nts]
plus_pos = list(np.where(np.in1d(chrom_seq, nts) == True)[0])
minus_pos = list(np.where(np.in1d(chrom_seq, rev_nts) == True)[0])
test_d['+'][chrom] = plus_pos
test_d['-'][chrom] = minus_pos
e = time.time()
print 'The numpy version took {} minutes...'.format((e-s)/60)
</code></pre>
<p>酵母产量:</p>
<pre><code>The numpy version took 0.0309354345004 minutes...
</code></pre>
<p>人类的输出:</p>
<pre><code>The numpy version took 9.86174853643 minutes...
</code></pre>
<p>为什么numpy版本在更大的人类基因组中失去了它的性能优势?有没有更快的方法?我尝试使用<code>multiprocessing.Pool</code>实现一个版本,但这比任何一个版本都慢:</p>
<pre><code>def getTestHelper(chrom_seq, nts, rev_comp):
rev_comps = { 'A' : 'T', 'T' : 'A', 'G' : 'C', 'C' : 'G', 'N' : 'N'}
if rev_comp:
nts = [rev_comps[nt] for nt in nts]
else:
nts = list(nts)
chrom_seq = np.array(list(chrom_seq))
mask = np.in1d(chrom_seq, nts)
pos_l = list(np.where(mask == True)[0])
return pos_l
s = time.time()
pool = Pool(4)
plus_pos = pool.map(functools.partial(getTestHelper, nts=nts, rev_comp=False), seq_d.values())
minus_pos = pool.map(functools.partial(getTestHelper, nts=nts, rev_comp=True), seq_d.values())
e = time.time()
print 'The parallel version took {} minutes...'.format((e-s)/60)
</code></pre>
<p>我还没有在人类基因组上运行过,但是酵母版本的速度较慢:</p>
<pre><code>The parallel version took 0.652778700987 minutes...
</code></pre>