игра хаос для последовательностей ДНК - PullRequest
Новая
       53

игра хаос для последовательностей ДНК

7 голосов
/ 04 ноября 2011

Я попробовал код mathematica для создания игры хаоса для последовательностей ДНК, размещенных по этому адресу: http://facstaff.unca.edu/mcmcclur/blog/GeneCGR.html

что-то вроде этого: 

genome = Import["c:\data\sequence.fasta", "Sequence"];
genome = StringReplace[ToString[genome], {"{" -> "", "}" -> ""}];
chars = StringCases[genome, "G" | "C" | "T" | "A"];
f[x_, "A"] := x/2;
f[x_, "T"] := x/2 + {1/2, 0};
f[x_, "G"] := x/2 + {1/2, 1/2};
f[x_, "C"] := x/2 + {0, 1/2};
pts = FoldList[f, {0.5, 0.5}, chars];
Graphics[{PointSize[Tiny], Point[pts]}]

последовательность фаста, которая у меня есть, это просто последовательность букв, как AACCTTTGATCAAA и генерируемый граф выглядит так:

enter image description here

код отлично работает с небольшими последовательностями, но когда я хочу поместить огромную последовательность, например, почти 40 Мб хромосомы, программа занимает много времени и отображает только черный квадрат, поэтому анализ невозможно. Можно ли улучшить вышеупомянутый код, чтобы квадрат, в котором он отображался, был бы больше? Кстати, квадрат должен быть только квадратной единицей. Заранее спасибо за помощь

1 Ответ

12 голосов
/ 04 ноября 2011

Сводная информация о приращениях, приведенных ниже:

Это значительно ускорит вычисление координат точек с использованием скомпилированного кода (в 50 раз, исключая вычисления shifts): 

shifts = chars /. {"A" -> {0., 0.}, "T" -> {.5, 0.}, "G" -> {.5, .5}, "C" -> {0, .5}};
fun1d = Compile[{{a, _Real, 1}}, FoldList[#/2 + #2 &, .5, a], CompilationTarget -> "C"]
pts = Transpose[fun1d /@ Transpose[shifts]];

Узким местом в вашем коде является рендеринг графики, мы вместо того, чтобы строить каждую точку, мы визуализируем плотность точек: 

threshold = 1;
With[{size = 300}, 
 Image[1 - UnitStep[BinCounts[pts, 1/size, 1/size] - threshold]]
]

Область будет окрашена в черный цвет, если онаимеет не менее threshold баллов.size - размер изображения.Выбрав большой размер или большой порог, вы можете избежать «проблемы черного квадрата».

Мой оригинальный ответ с более подробной информацией:

На моемдовольно устаревшая машина, код не очень медленный. 

chars = RandomChoice[{"A", "T", "C", "G"}, 800000];

f[x_, "A"] := x/2;
f[x_, "T"] := x/2 + {1/2, 0};
f[x_, "G"] := x/2 + {1/2, 1/2};
f[x_, "C"] := x/2 + {0, 1/2};
Timing[pts = FoldList[f, {0.5, 0.5}, chars];]
Graphics[{PointSize[Tiny], Point[pts]}]

Я получаю время 6,8 секунды, которое можно использовать, если вам не нужно запускать его много раз в цикле (если это не достаточно быстро дляВаш вариант использования и устройство, пожалуйста, добавьте комментарий, и мы постараемся ускорить его.)

К сожалению, рендеринг графики занимает намного больше времени (36 секунд), и я не знаю, есть ливсе, что вы можете с этим поделать.Отключение сглаживания может немного помочь , в зависимости от вашей платформы, но не сильно: Style[Graphics[{PointSize[Tiny], Point[pts]}], Antialiasing -> False] (для меня это не так).Это давнее раздражение для многих из нас.

Что касается всей графики, являющейся черной, вы можете изменить ее размер с помощью мыши и сделать ее больше.В следующий раз, когда вы оцените свое выражение, выходной рисунок запомнит его размер.Или просто используйте ImageSize -> 800 в качестве опции Graphics.Учитывая плотность пикселей на экранах, единственное другое решение, которое я могу придумать (которое не включает изменение размера графики), состоит в том, чтобы представить плотность пикселей с использованием оттенков серого и построить плотность.

РЕДАКТИРОВАТЬ:

Вот как вы можете построить плотность (это также гораздо быстрее вычислить и визуализировать, чем точечный график!): 

With[{resolution = 0.01}, 
 ArrayPlot@BinCounts[pts, resolution, resolution]
]

Играть сразрешение, чтобы сделать сюжет хорошим.

Для моего примера с произвольной последовательностью это дает только серый график.Для ваших данных генома это, вероятно, даст более интересную схему.

РЕДАКТИРОВАТЬ 2:

Вот простой способ ускорить функцию с помощью компиляции:

Сначала замените символы векторами сдвига (это нужно сделать только один раз для набора данных, затем вы можете сохранить результат): 

arr = chars /. {"A" -> {0., 0.}, "T" -> {.5, 0.}, "G" -> {.5, .5}, "C" -> {0, .5}};

Затем скомпилируем нашу функцию: 

fun = Compile[{{a, _Real, 2}}, FoldList[#/2 + #2 &, {.5, .5}, a], 
 CompilationTarget -> "C"]

Удалите CompilationTarget, если ваша версия Mathematica более ранняя, чем 8, или у вас не установлен компилятор C. 

fun[arr]; // Timing

дает мне 0,6 секунды, что является мгновенным 10-кратным ускорением.

РЕДАКТИРОВАТЬ 3:

Возможно еще 5-кратное ускорение по сравнению с вышеупомянутой скомпилированной версией, благодаря отсутствию некоторых обратных вызовов ядра в скомпилированной функции (я проверял вывод компиляции, используя CompilePrint чтобы придумать эту версию - иначе это не очевидно почему это быстрее): 

fun1d = Compile[{{a, _Real, 1}}, FoldList[#/2 + #2 &, .5, a], 
  CompilationTarget -> "C"]

arrt = Transpose[arr];
Timing[result = fun1d /@ arrt;]
pts = Transpose[result];

На моем компьютере это выполняется за 0,11 секунды.На более современной машине он должен завершиться через несколько секунд даже для набора данных 40 МБ.

Я разделил транспонирование на отдельные входы, потому что в этот момент время выполнения fun1d начинает сравниваться свремя работы Transpose.

...