Я хочу преобразовать приведенный ниже код в список, чтобы оптимизировать его для скорости. Он использует два вложенных цикла, а входные Motifs - список строк, содержащих A C G и T.
Я попытался заменить вложенные для циклов внизу функции для этого понимания списка:
Это не работает, возможно потому, что мне не разрешено присваивать значение + = 1 в понимании списка. Как я могу обойти это?
Ответ 3
Я предлагаю метод numpy, предоставленный @Kasramvs, если производительность скорости действительно имеет значение.
Кроме того, подсчет символов не подходит даже для современных компьютеров, и, возможно, вы можете играть с некоторыми трюками об индексировании/хэшировании входных строк перед подсчетом. Например, поскольку каждая строка имеет только 4 символа, и каждый char содержит только 4 возможных буквы: "A", "C", "G", "T", поэтому он может легко представлять каждый из всех "ACGT 'комбинации от" AAAA "до" TTTT" с числом, хешем или таинственным кодом. Объем комбинаций должен быть равен или меньше 4x4x4x4 = 256 различных номеров здесь.
И затем пересчитайте код. Например, каждый раз, когда вы видите "AAAA", а затем считаете его 0x0 в списке python или массивом numpy, см. "AAAC" и считайте 0x1 и наоборот. После этого вы получите массив binning с индексами от 0x0 ~ 0xFF (255) и с соответствующими вхождениями, правильно? Теперь помните, что один 0x0 означает A: {1, 1, 1, 1} в вашем случае или семь 0x1 для A: {7, 7, 7, 0} вместе с C: {0, 0, 0, 7}... Суммируйте их все тщательно и это результат.
Эти трюки, возможно, в этом случае могут значительно повысить скорость работы. Скорость выигрывает в двух факторах: первая заключается в том, что теперь компьютеры имеют дело с числами вместо символов, в то время как числа намного проще сортировать, считать, группировать, разбивать или индексировать, чем символы; второй - из-за гораздо более высокой скорости попадания в кэш по характеру, поскольку трассировка памяти в этих трюках значительно сокращается.
Надеюсь, это поможет.:)
Ну, добавьте некоторые коды, как показано ниже, как четкую ссылку.
Прежде всего, импорт:
import itertools
import numpy as np
И нижеприведенная функция encode_sequence()
должна быть достаточно быстрой, если dict t02
не слишком велик, скажем, обычно менее 1M пар ключ-значение:
def encode_sequence(tsq):
t00 = itertools.product('ACGT', repeat=4)
t01 = np.array(list(t00)).view('|U4').ravel()
t02 = dict(zip(t01, np.arange(len(t01))))
t03 = [t02[i] for i in tsq]
return t03
И используйте следующий фрагмент, чтобы сгенерировать тензор map_decode
, чтобы представить материал о `` counting code ''... Кроме того, есть математический трюк, называемый расширенным преобразованием матрицы под этой частью, говоря, что он преобразует ll0
и 'ACGT'
таким же образом, чтобы охватить map_decode
для последующего использования.
ll0 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
map_decode = np.array(list(itertools.product(ll0, repeat=4)))
Подайте тестовую последовательность и переведите,
test_seq = ('ACGG', 'CAGT', 'CCGA', 'AAAT', 'TTGC', 'GCAT', 'ACGG', 'AAAT')
u01 = encode_sequence(test_seq)
Подсчитайте вхождения; что блок ниже должен быть основным источником усиления скорости, потому что компьютер хорошо справляется с числами в u01
,
p01, q01 = np.unique(u01, return_counts=True)
В конце концов, сгенерируйте вывод... Здесь немного сложно, например p01
является отсортированным хеш-кодом test_seq
и q01
, действительно, соответствующими подсчетами, а map_decode
служит тем, что я сказал, что тензор отображает хэш-код в p01
на другой вектор, который мы хотим, например, отображение 0x0 (или 'AAAA') на A: [1, 1, 1, 1], C: [0, 0, 0, 0], G: [0, 0, 0, 0] и T: [0, 0, 0, 0]. Отображаемый map_decode[p01]
таким образом взвешивается счетчиками q01
и готов к сумме для отчета:
np.sum(map_decode[p01]*q01[:, None, None], axis=0).T
И он говорит:
array([[4, 3, 3, 1],
[2, 4, 0, 1],
[1, 0, 5, 2],
[1, 1, 0, 4]])
что эквивалентно A: {4, 3, 3, 1}, C: {2, 4, 0, 1}, G: {1, 0, 5, 2} и T: {1, 1, 0, 4}. Проверьте, соответствует ли он ответу.
Это реализация numpy; он не содержит явного цикла в основном теле. И более того, encode_sequence()
может быть заменен некоторыми оффлайн-подготовками входов для повышения производительности заблаговременно. Хотя у меня не было показателя скорости вышеприведенных фрагментов, я думаю, что их нужно ускорить до определенной степени.:)
Хорошо, обсудим, что произойдет, если есть длинные строки.
Мы используем эту последовательность в качестве примера,
test_seq0 = ((
'A'*40, 'A'*40, 'A'*40, 'C'*40, 'C'*40,
'C'*40, 'C'*40, 'G'*40, 'G'*40, 'T'*40
))*4
test_seq0
содержит 40 строк, каждая строка содержит 40 символов. Это выглядит так:
In: len(test_seq0), test_seq0
Out: (40,
('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
... skip 30 lines ...
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC',
'GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG',
'GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG',
'TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT'))
Довольно забавный вид, верно?
Затем мы должны установить encode_sequence()
для длинной строки,
def encode_sequence_longstring(tsq_np):
t00 = itertools.product('ACGT', repeat=4)
t01 = np.array(list(t00)).view('|U4').ravel()
t02 = dict(zip(t01, np.arange(len(t01))))
t03 = np.empty_like(tsq_np, dtype=np.uint)
t03.ravel()[:] = [t02[i] for i in tsq_np.ravel()]
return t03
Будьте осторожны, что tsq_np
здесь больше не является простым списком строк. Постфикс _np
означает теперь массив numpy.
И разделите оригинал test_seq0
на несколько шагов,
In: v01 = np.asarray(test_seq0).view('|U4').reshape(-1, int(40/4))
In: v01
Out:
array([['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA', 'AAAA'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
... skip 30 lines ...
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC', 'CCCC'],
['GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG'],
['GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG', 'GGGG'],
['TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT', 'TTTT']],
dtype='<U4')
Еще один забавный вид в v01
.:)
И используйте v01
для вычисления хэш-кодов u02
следующим образом. Это связано с некоторыми соглашениями о множественных числах вокруг этих переменных и функций. Просто привыкнуть к этим причудливым трюкам; они того стоят,
In: u02 = encode_sequence_longstring(v01)
In: u02
Out:
array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
... skip 30 lines ...
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[ 85, 85, 85, 85, 85, 85, 85, 85, 85, 85],
[170, 170, 170, 170, 170, 170, 170, 170, 170, 170],
[170, 170, 170, 170, 170, 170, 170, 170, 170, 170],
[255, 255, 255, 255, 255, 255, 255, 255, 255, 255]],
dtype=uint64)
По наблюдениям вы можете сказать, что u02
на самом деле является отображением 1-to-1 для v01
. Он просто отображает каждый "AAAA" в 0x0, как ожидалось.
Теперь карта u02
содержит всю необходимую информацию относительно test_seq0
. Извлеките его из u02
с помощью numpy,
s01 = np.empty((4, 0))
for u03 in u02.T:
p02, q02 = np.unique(u03, return_counts=True)
s02 = np.sum(map_decode[p02]*q02[:, None, None], axis=0).T
s01 = np.hstack((s01, s02))
Там правило большого пальца для собственного цикла python: вы можете использовать его, но использовать его вне любой чувствительной к производительности области. Но это требует опыта, чтобы судить о ситуации в руке.
Теперь s01
- ожидаемый отчет следующим образом:
In: s01
Out:
array([[ 12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.,
12., 12., 12., 12., 12., 12., 12.],
[ 16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16., 16., 16., 16., 16.,
16., 16., 16., 16., 16., 16., 16.],
[ 8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8., 8., 8., 8., 8.,
8., 8., 8., 8., 8., 8., 8.],
[ 4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4., 4., 4., 4., 4.,
4., 4., 4., 4., 4., 4., 4.]])
Прочитайте 4 строки сверху вниз, и они - "A", "C", "G", "T" соответственно.
В то же время я пробовал 40x10x4000 test_seq0
, как это,
test_seq0 = ((
'A'*40, 'A'*40, 'A'*40, 'C'*40, 'C'*40,
'C'*40, 'C'*40, 'G'*40, 'G'*40, 'T'*40
))*4000
В отчете говорится:
array([[ 12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000., 12000., 12000.,
12000., 12000., 12000., 12000., 12000.],
[ 16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000., 16000., 16000.,
16000., 16000., 16000., 16000., 16000.],
[ 8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000., 8000., 8000.,
8000., 8000., 8000., 8000., 8000.],
[ 4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000., 4000., 4000.,
4000., 4000., 4000., 4000., 4000.]])
и он заканчивается менее чем за 1 секунду (нажмите ENTER и сделайте это) на моем MacBook Pro, который не является монстром.:)