matlab - 如何绘制 DNA 序列的基因图，例如 ATGCCGCTGCGC？

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），该文章还给出了 IFS 的混沌游戏算法，使用 ( DNA 序列的四个碱基。可以在此处找到详细说明（原始文章）。

该方法可用于生成如下图。只是为了它，我已经包括（在 RHS 面板中）用相应的互补 DNA 链 (cDNA) 生成的图。

小鼠线粒体 DNA (LHS) 及其 互补链 (cDNA) (RHS)。

在此处输入图像描述

这些图是从 GenBank 标识符 gi|342520 生成的。该序列包含 16295 个碱基。

（Jose Manuel Gutiérrez 使用的示例之一。如果有人感兴趣，可以从 gi|1262342 生成人类等效图）。

人 β 珠蛋白区域(LHS) 及其 cDNA (RHS)

在此处输入图像描述

生成自 gi|455025| （该示例使用了我的 Mark McClure）。该序列包含 73308 个碱基

有相当有趣的情节！这种图的（有时）分形性质是已知的，但在 LHS 与 RHS（cDNA）版本中明显的对称性非常令人惊讶（至少对我而言）。

好消息是，通过直接导入序列（例如来自 Genbank），然后使用 Mma 的强大功能，可以很容易地生成任何 DNA 序列的此类图。
所有你需要它的入藏号！（“未知”核苷酸，如“R”可能需要被删除）（我使用的是 Mma v7）。

The Original Implimenation (略有修改) (by Jose Manuel Gutiérrez)

重要更新

在 Mark McClure 的建议下，我已更改Point/@Orbit[s, Union[s]]为Point@Orbit[s, Union[s]].

这大大加快了速度。请参阅下面的 Mark 评论。

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

第二组图以类似的方式从 gi|455025 生成。请注意，序列很长！

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

最后一个例子（包含 265922 bp），也显示出迷人的“分形”对称性。（这些是用AbsolutePointSize[1]in生成的IFSPlot）。

fasta 文件的第一行：

{>gi|328530803|gb|AFBL01000008.1| 放线菌属口腔分类群 170 str。F0386 A_spOraltaxon170F0386-1.0_Cont9.1，全基因组鸟枪法序列}

在此处输入图像描述

相应的 cDNA 图在 RHS 上再次以蓝色显示

最后，Mark 的方法还给出了非常漂亮的绘图（例如 gi|328530803），并且可以作为笔记本下载。

score 28 · Accepted Answer

并不是说我真的理解你想要的“图表”，但这是一种字面解释。

以下代码都不一定是最终形式。在我尝试改进任何东西之前，我想知道这是否正确。

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，或者在 1990 年由 Joel Jefferey 撰写的论文“基因结构的混沌游戏表示”中描述的基因序列的所谓混沌游戏表示。这是 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]]

以下是如何将其应用于从 Mathematica 的GenomeData命令中获取的基因序列：

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

score 0 · Accepted Answer

这个问题似乎已经得到了很好的回答，但我想我要补充一点，你所描述的内容之前已经在DNA 序列的数字表示方法的集合中以DNA的名义发表，这些方法在我们的预印本中进行了讨论。

事实证明，DNA 行走在实践中并不是很有用，但可以进行直观的可视化。我手头没有，但我想我的同事会非常乐意提供用于生成下图的 Matlab 代码。

matlab - 如何绘制 DNA 序列的基因图，例如 ATGCCGCTGCGC？

6 回答 6

3D 选项

Related

Reference