7

我已经尝试过使用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 的染色体，该程序需要很多时间并且只显示一个黑色方块，因此无法分析。是否可以改进上述代码，使其显示的正方形更大？顺便说一下，正方形必须只是正方形单位。提前感谢您的帮助

4

1 回答 1

12

以下增量编辑的摘要：

这将使您在使用编译代码（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。

于 2011-11-04T13:04:58.400 回答