wolfram-mathematica – DNA序列的混沌游戏

我已经尝试了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的染色体时,该程序需要花费大量时间并且只显示黑色方块,因此无法进行分析.
是否有可能改进上述代码,以便显示它的方格会更大?,方式必须只是方形单位.
感谢您的帮助

以下增量编辑摘要:

通过使用编译代码(50x不包括计算班次),这将为您提供相当快的计算点坐标的速度:

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]]
]

如果区域至少具有阈值点,则该区域将被涂成黑色.尺寸是图像尺寸.通过选择大尺寸或大阈值,您可以避免“黑方问题”.

我的原始答案更详细:

在我相当过时的机器上,代码不是很慢.

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秒),我不知道你能做些什么.禁用抗锯齿可能会有所帮助,具体取决于您的平台,但不多:样式[Graphics [{PointSize [Tiny],Point [pts]}],Antialiasing – >错](对我来说没有).这对我们许多人来说是一个长期的烦恼.

关于整个图形为黑色,您可以使用鼠标调整大小并使其更大.下次评估表达式时,输出图形将记住其大小.或者只使用ImageSize – > 800作为图形选项.考虑到屏幕的像素密度,我能想到的唯一其他解决方案(不涉及调整图形大小)将使用灰度阴影表示像素密度,并绘制密度.

编辑:

这就是你可以绘制密度的方法(计算和渲染的速度比点图快得多!):

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

相关文章
相关标签/搜索