matlab - ATGCCGCTGCGCと言うDNA配列の遺伝子グラフをプロットする方法は？

Question

ウイルスの塩基対配列が2k塩基対であるため、ウイルスのDNA配列に基づいてランダムウォークを生成する必要があります。シーケンスは「ATGCGTCGTAACGT」のようになります。パスは、Aの場合は右に、Tの場合は左に、Gの場合は上に、Cの場合は下に曲がる必要があります。この目的でMatlab、Mathematica、またはSPSSを使用するにはどうすればよいですか？

score 33 · Accepted Answer

遺伝子配列のカオスゲーム表現に関するMark McClure のブログを以前は知りませんでしたが、Jose Manuel Gutiérrez の記事 (The Mathematica Journal Vol 9 Issue 2) を思い出しました。この記事では、( DNA 配列の 4 つの塩基。詳細な説明はこちら(元の記事) にあります。

このメソッドを使用して、次のようなプロットを作成できます。念のために、対応する相補的 DNA 鎖 (cDNA) で生成されたプロットを (RHS パネルに) 含めました。

マウスミトコンドリア DNA (LHS) とその 相補鎖 (cDNA) (RHS)。

ここに画像の説明を入力

これらのプロットは、GenBank Identifier gi|342520 から生成されました。この配列には 16295 塩基が含まれています。

(Jose Manuel Gutiérrez によって使用された例の 1 つです。興味のある方は、gi|1262342 から人間に相当するプロットを生成できます)。

ヒト β グロビン領域(LHS) とその cDNA (RHS)

ここに画像の説明を入力

gi|455025| から生成 (この例では、私の Mark McClure を使用しました)。シーケンスには 73308 塩基が含まれます

かなり興味深いプロットがあります！このようなプロットの (場合によっては) フラクタルな性質は知られていますが、LHS と RHS (cDNA) のバージョンで明らかな対称性は非常に驚くべきものでした (少なくとも私にとっては)。

良い点は、任意の DNA 配列のこのようなプロットは、シーケンスを (たとえば、Genbank から)直接インポートし、Mma の機能を使用することによって非常に簡単に生成できることです。
必要なのはアクセッション番号だけです! (「R」などの「未知の」ヌクレオチドはザッピングする必要がある場合があります) (私は Mma v7 を使用しています)。

The Original Implimenation (わずかに修正) (Jose Manuel Gutiérrez 作)

重要な更新

Mark McClure の助言により、私はに変更Point/@Orbit[s, Union[s]]しましたPoint@Orbit[s, Union[s]]。

これにより、処理が大幅に高速化されます。以下のマークのコメントを参照してください。

Orbit[s_List, {a_, b_, c_, d_}] := 
  OrbitMap[s /. {a -> {0, 0}, b -> {0, 1}, c -> {1, 0}, 
     d -> {1, 1}}];
OrbitMap = 
  Compile[{{m, _Real, 2}}, FoldList[(#1 + #2)/2 &, {0, 0}, m]];
IFSPlot[s_List] := 
 Show[Graphics[{Hue[{2/3, 1, 1, .5}], AbsolutePointSize[2.5], 
    Point @ Orbit[s, Union[s]]}], AspectRatio -> Automatic, 
  PlotRange -> {{0, 1}, {0, 1}}, 
  GridLines -> {Range[0, 1, 1/2^3], Range[0, 1, 1/2^3]}]

これにより、青いプロットが得られます。緑の場合、Hue[] を Hue[{1/3,1,1,.5}] に変更します。

次のコードは、最初のプロットを生成します (マウスのミトコンドリア DNA 用)。

 IFSPlot[Flatten@
      Characters@
       Rest@Import[
         "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=\
    nucleotide&id=342520&rettype=fasta&retmode=text", "Data"]]

cDNAプロットを取得するために、次の変換ルールを使用しました（また、色相設定も変更しました）

IFSPlot[    ....   "Data"] /. {"A" -> "T", "T" -> "A", "G" -> "C", 
   "C" -> "G"}]

Sjoerd C. de Vriesとtelefunkenvf14のおかげで、NCBI サイトから配列を直接インポートすることができました。

わかりやすくするために、物事を少し分割します。

シーケンスをインポートする

mouseMitoFasta=Import["http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=342520&rettype=fasta&retmode=text","Data"];

Mathematica J.の元の記事でシーケンスをインポートする方法は古いものです。

素敵なチェック

はじめ@mouseMitoFasta

出力：

{>gi|342520|gb|J01420.1|MUSMTCG Mouse mitochondrion, complete genome}

塩基リストの生成

mouseMitoBases=Flatten@Characters@Rest@mouseMitoFasta

さらにいくつかのチェック

{Length@mouseMitoBases, Union@mouseMitoBases,Tally@mouseMitoBases}

出力：

{16295,{A,C,G,T},{{G,2011},{T,4680},{A,5628},{C,3976}}}

2 番目のプロットセットは、gi|455025 から同様の方法で生成されました。シーケンスが長いことに注意してください。

{73308,{A,C,G,T},{{G,14785},{A,22068},{T,22309},{C,14146}}}

最後の例(265922 bp を含む) も、魅力的な「フラクタル」対称性を示しています。(これらは in で生成されAbsolutePointSize[1]ましたIFSPlot)。

fasta ファイルの最初の行:

{>gi|328530803|gb|AFBL01000008.1| 放線菌sp。口腔分類群 170 str. F0386 A_spOraltaxon170F0386-1.0_Cont9.1, 全ゲノムショットガンシーケンス}

ここに画像の説明を入力

対応する cDNA プロットは、RHS に青色で表示されます。

最後に、Mark の方法でも非常に美しいプロットが得られ (たとえば gi|328530803 を使用)、ノートブックとしてダウンロードできます。

score 28 · Accepted Answer

私はあなたが望む「グラフ」を本当に理解しているわけではありませんが、ここに1つの文字通りの解釈があります。

次のコードはいずれも、必ずしも最終的な形式ではありません。何かを洗練しようとする直前に、これが正しいかどうか知りたいです。

rls = {"A" -> {1, 0}, "T" -> {-1, 0}, "G" -> {0, 1}, "C" -> {0, -1}};
Prepend[Characters@"ATGCGTCGTAACGT" /. rls, {0, 0}];
Graphics[Arrow /@ Partition[Accumulate@%, 2, 1]]

ここに画像の説明を入力してください

Prepend[Characters@"TCGAGTCGTGCTCA" /. rls, {0, 0}];
Graphics[Arrow /@ Partition[Accumulate@%, 2, 1]]

ここに画像の説明を入力してください

3Dオプション

i = 0;
Prepend[Characters@"ATGCGTCGTAACGT" /. rls, {0, 0}];
Graphics[{Hue[i++/Length@%], Arrow@#} & /@ 
  Partition[Accumulate@%, 2, 1]]

ここに画像の説明を入力してください

i = 0;
Prepend[Characters@"ATGCGTCGTAACGT" /. 
    rls /. {x_, y_} :> {x, y, 0.3}, {0, 0, 0}];
Graphics3D[{Hue[i++/Length@%], Arrow@#} & /@ 
  Partition[Accumulate@%, 2, 1]]

ここに画像の説明を入力してください

必要なものがわかったので、最初の関数のパッケージバージョンを次に示します。

genePlot[s_String] :=
 Module[{rls},
  rls =
   {"A" -> { 1, 0},
    "T" -> {-1, 0},
    "G" -> {0,  1},
    "C" -> {0, -1}};
  Graphics[Arrow /@ Partition[#, 2, 1]] & @
   Accumulate @ Prepend[Characters[s] /. rls, {0, 0}]
]

次のように使用します。

genePlot["ATGCGTCGTAACGT"]

score 16 · Accepted Answer

CGR、またはJoel Jefferey による1990 年の論文「遺伝子構造のカオスゲーム表現」で説明されている遺伝子配列のいわゆるカオスゲーム表現について話しているように思えます。Mathematica での実装は次のとおりです。

cgrPic[s_String] := Module[
  {},
  chars = StringCases[s, "G"|"A"|"T"|"C"];
  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];
  ListPlot[pts, AspectRatio -> Automatic]]

GenomeDataMathematica のコマンドから取得した遺伝子配列に適用する方法は次のとおりです。

cgrPic[GenomeData["FAT4", "FullSequence"]]

score 6 · Accepted Answer

シーケンスがすでに*)整数配列にマッピングされていると仮定するとS、実際の動きの計算はルールに基づいて簡単ですR:

R =
   1  -1   0   0
   0   0   1  -1
S =
   1   2   3   4   3   2   4   3   2   1   1   4   3   2
T= cumsum(R(:, S), 2)
T =
   1   0   0   0   0  -1  -1  -1  -2  -1   0   0   0  -1
   0   0   1   0   1   1   0   1   1   1   1   0   1   1

*) 実際のシーケンスについては、さらに詳しく説明する必要があります。それは単一の文字列として表されているのでしょうか、それともセル配列として表されているのでしょうか?
編集:
シーケンスが文字列として表されていると仮定すると、次のSような整数シーケンスにマップします。

r= zeros(1, 84);
r(double("ATGC"))= [1 2 3 4];
S= r(double("ATGCGTCGTAACGT"))

そしてそれをプロットするには：

plot([0 T(1, :)], [0 T(2, :)], linespec)

はlinespec目的の回線仕様です。

score 6 · Accepted Answer

このようなものを試すこともできます...

RandomDNAWalk[seq_, path_] := 
 RandomDNAWalk[StringDrop[seq, 1], 
  Join[path, getNextTurn[StringTake[seq, 1]]]];

RandomDNAWalk["", path_] := Accumulate[path];

getNextTurn["A"] := {{1, 0}};
getNextTurn["T"] := {{-1, 0}};
getNextTurn["G"] := {{0, 1}};
getNextTurn["C"] := {{0, -1}};

ListLinePlot[
 RandomDNAWalk[
  StringJoin[RandomChoice[{"A", "T", "C", "G"}, 2000]], {{0, 0}}]]

matlab - ATGCCGCTGCGCと言うDNA配列の遺伝子グラフをプロットする方法は？

6 に答える 6

3Dオプション

Related

Reference