7

このアドレスに投稿された DNA シーケンスのカオス ゲームを作成するための 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]}]

私が持っている fasta シーケンスは AACCTTTGATCAAA のような文字のシーケンスであり、生成されるグラフは次のようになります。

ここに画像の説明を入力

コードは小さなシーケンスでは問題なく動作しますが、たとえば約 40Mb の染色体など、巨大なシーケンスを配置したい場合、プログラムは多くの時間を要し、黒い四角が表示されるだけなので分析できません。前述のコードを改善して、表示される正方形が大きくなるようにすることはできますか?ところで、正方形は正方形単位のみである必要があります。事前にご協力いただきありがとうございます

4

1 に答える 1

12

以下の増分編集の概要:

これにより、コンパイルされたコードを使用してポイント座標を計算する際の速度が大幅に向上します ( 計算を除くと 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：

コンパイルを使用して関数を高速化する簡単な方法を次に示します。

まず、文字をシフト ベクトルに置き換えます (データセットに対して 1 回だけ実行する必要があり、結果を保存できます)。

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

CompilationTargetMathematica のバージョンが 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 に答える