DNA序列的混沌游戏

Question

我已经尝试了mathematica代码,用于在这个地址中发布DNA序列的混沌游戏: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]}]

我所拥有的fasta序列只是像AACCTTTGATCAAA这样的字母序列,要生成的图形如下:

代码适用于小序列,但是当我想要放置一个巨大的序列,例如几乎40Mb的染色体时,该程序需要花费大量时间并且只显示黑色方块,因此无法进行分析.是否有可能改进上述代码,以便显示它的方格会更大？,方式必须只是方形单位.感谢您的帮助

Answer 1

增量编辑摘要如下:

通过使用编译代码(50x不包括计算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 MB的数据集,它也应该在几秒钟内完成.

我将换位拆分为单独的输入,因为此时fun1d开始的运行时间与运行时间相当Transpose.