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

以下增量修改摘要：

这将使您在使用编译后的代码计算点坐标时显着加快（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"]

如果您的 Mathematica 版本早于 8 或您没有安装 C 编译器，请删除 CompilationTarget。

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 的运行时间相当。