如何绘制DNA序列的基因图表说ATGCCGCTGCGC？
imb*_*ene 34 matlab wolfram-mathematica bioinformatics dna-sequence
我需要根据病毒的DNA序列生成随机游走,给出其2k碱基对的碱基对序列.序列看起来像"ATGCGTCGTAACGT".路径应向右转为A,向左转为T,向上转为G,向下转向C.如何为此目的使用Matlab,Mathematica或SPSS？
tom*_*omd 30
我以前不知道马克麦克卢尔关于基因序列的混沌游戏表现的博客,但它让我想起了JoseManuelGutiérrez撰写的一篇文章(The Mathematica Journal Vol 9 Issue 2),它也为IFS提供了一个混沌游戏算法( DNA序列的四个碱基.可在此处找到详细描述(原始文章).
该方法可用于产生如下的图.只是为了它的地狱,我已经包括(在RHS面板中)用相应的互补DNA链(cDNA)产生的图.
- 小鼠线粒体DNA(LHS)及其 互补链(cDNA)(RHS).
这些图是从GenBank Identifier gi | 342520生成的.该序列含有16295个碱基.
(JoseManuelGutiérrez使用的一个例子.如果有人感兴趣,可以从gi | 1262342生成人类等效物的图.
- 人β珠蛋白区(LHS)及其cDNA(RHS)
源自gi | 455025 | (这个例子使用了我的Mark McClure).该序列含有73308个碱基
有很有趣的情节!这种图的(有时)分形性质是已知的,但LHS与RHS(cDNA)版本中明显的对称性是非常令人惊讶的(至少对我而言).
好的是,通过直接输入序列(来自Genbank),然后使用Mma的功能,可以非常容易地生成任何DNA序列的这种图 .
您需要的所有入藏号码!('未知'核苷酸如"R"可能需要被破坏)(我正在使用Mma v7).
The Original Implimenation(略有修改) (作者JoseManuelGutiérrez)
重要更新
在马克·麦克卢尔的提醒,我已经改变了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图,我使用了以下转换规则(并且还更改了Hue设置)
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),可以下载为笔记本.
- 谢谢你指出那篇文章.了解如何从该站点自动导入特别好.我有两条评论:首先,我提到的网站似乎可供我使用.您确定URL中没有拼写错误吗？您可以通过从页面顶部的链接下载该帖子的笔记本版本来避免该潜在问题.其次,通过将"Point/@ Orbit [s,Union [s]]"更改为"Point @ Orbit [s,Union [s]]",可以大大加快Jose的代码速度.在该文章出现之后,这是在V6中引入的Point的一种新的重要语法. (2认同)
Mr.*_*ard 28
并不是说我真的理解你想要的"图形",但这是一个字面解释.
以下代码中的任何一个都不一定是最终形式.在我尝试改进任何东西之前,我想知道这是否正确.
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"]
Mar*_*ure 16
听起来你在谈论CGR,或者所谓的基因序列的混沌游戏表示.几个月前我在博客上写了这篇文章:http: //facstaff.unca.edu/mcmcclur/blog/GeneCGR.html
假设序列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所需的线规格在哪里.
你也可以试试这样的......
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}}]]