转自:http://blog.sina.com.cn/s/blog_5ecfd9d90100cm03.html~type=v5_one&label=rela_prevarticle
matlab工具箱中的生物信息学工具箱(bioinformatics toolbox)功能还算全面,虽然没有什么突出的优点。包含了:蛋白和核酸分析,系统发育分析以及基因芯片分析等功能。
1、序列分析:以下,带有>>的为命令行!!
例如,我们要查询 RefSeq的一条序列NM_000799,
>> a = getgenbank('NM_000799','ToFile','testbank'); %利用getgenbank函数,可以设置所需信息的各种参数。
%getgenbank('accession',TOFILE',FILENAME):可以将信息写入文件;
%getgenbank('accession','SEQUENCEONLY',true);可以只提取序列信息;
%getgenbank('accession','PARTIALSEQ',SEQPARAMS);可以只提取一部分序列,SEQPARAMS为[N,M];
%getgenbank('accession','fileformat',fmt);fmt可以为genbank或者fasta。
>>ntseq = a.Sequence %这个是结构数组,就是在数组变量名后面加一个圆点,例如array.substructure
>>ntdensity(ntseq) %图示核酸序列各种核苷酸含量;可以设置窗口大小
%ntdensity(ntseq,'window',windowvalue)
>>basecount(ntseq,'chart','bar') %计算核苷酸个数,也可以图示化basecount%(ntseq,'chart','value') value可以设为bar或者pie
>> basecount(seqrcomplement(ntseq))
%它的反义链的核苷酸个数
>> codoncount(ntseq) %计算密码子使用频率
%codoncount(SeqNT, ...'Frame', FrameValue, ...) FrameValue设置阅读框数值,1,2或者3
%codoncount(SeqNT, ...'Reverse', ReverseValue, ...)设置互补链,数值为true或者false
%codoncount(SeqNT, ...'Figure', FigureValue, ...)是否图示化,true或者false
%由于不知道阅读框(ORF)起始位置,所以有六种可能:
%写个循环:
>>for frame = 1:3 %每条链有三种可能
figure('color',[1 1 1]) %每次打开一个图形界面
subplot(2,1,1);%将图新划分成小块
codoncount(ntseq,'frame',frame,'figure',true);
title(sprintf('Codons for frame %d',frame));%添加标题
subplot(2,1,2);
codoncount(ntseq,'reverse',true,... %互补链
'frame',frame,...
'figure',true);
title(sprintf('Codons for reverse frame %d',frame)); %sprintf设置输出格式,%d十进制输出
end
>>f = seqshoworfs(ntseq)
%显示阅读框 seqshoworfs(SeqNT, ...'Frames', FramesValue, ...)阅读框设置,1,2或者3
%seqshoworfs(SeqNT, ...'GeneticCode', GeneticCodeValue, ...)不同的生物体密码子有所不同
>>StartIndex = find(f(1).Start ==196) %我们选第一个阅读框196位起始密码子
>>ND2Stop = f(1).Stop(StartIndex) %终止密码子
>> ND2Seq =ntseq(196:ND2Stop); %截选这一段序列
>>codoncount (ND2Seq) %密码子频率
%下面看一下,翻译过程:
>>ND2AASeq = nt2aa(ND2Seq);%翻译,可以设置密码子表
%SeqAA = nt2aa(..., 'GeneticCode', GeneticCodeValue, ...) 设置数值不一样,结果也不一样
>>aacount(ND2AASeq, 'chart','bar') %图示氨基酸频率
>>atomiccomp(ND2AASeq); %原子组成
>>molweight (ND2AASeq); %分子量大小
系统发生分析 -开始太傻了!!
1、首先建立malab结构,将要分析的各物种的信息输入:
>>data = {'German_Neanderthal' 'AF011222';
'Russian_Neanderthal' 'AF254446';
'European_Human' 'X90314' ;
'Mountain_Gorilla_Rwanda' 'AF089820';
'Chimp_Troglodytes' 'AF176766';
};%建立结构
2、开始准备序列,以备下面分析:
>>for ind = 1:5 %设置循环,因为上面结构一共五个物种
seqs(ind).Header = data{ind,1};%一样设置为结构
seqs(ind).Sequence = getgenbank(data{ind,2},...
'sequenceonly', true);
end
3、这儿就要讲一下,如何进行序列比对了:
D = seqpdist(Seqs)%可以是细胞数组,或者是包含序列结构向量,或者是矩阵,每一行一条序列;
D = seqpdist(Seqs, ...'Method', MethodValue, ...)% p距离或者(Default is Jukes-Cantor)或者序列比对得分。
D = seqpdist(Seqs, ...'Indels', IndelsValue, ...)% 插入缺失突变是否算分
D = seqpdist(Seqs, ...'Optargs', OptargsValue, ...)
D = seqpdist(Seqs, ...'PairwiseAlignment', PairwiseAlignmentValue, ...)%全局序列比对,长度不一则为TRUE,否则为FALSE
D = seqpdist(Seqs, ...'JobManager', JobManagerValue, ...)%需要进行Parallel Computing Toolbox分析
D = seqpdist(Seqs, ...'WaitInQueue', WaitInQueueValue, ...)
D = seqpdist(Seqs, ...'SquareForm', SquareFormValue, ...)%是否把结果输出方阵
D = seqpdist(Seqs, ...'Alphabet', AlphabetValue, ...)%'NT'标示序列为核酸 or 蛋白'AA'
D = seqpdist(Seqs, ...'ScoringMatrix', ScoringMatrixValue, ...)
% 算分矩阵
*
'PAM40'
*
'PAM250'
*
'DAYHOFF'
*
'GONNET'
*
'BLOSUM30' increasing by 5 up to 'BLOSUM90'
*
'BLOSUM62'
*
'BLOSUM100'
默认Default is:
*
'NUC44' (when AlphabetValue equals 'NT')
*
'BLOSUM50' (when AlphabetValue equals 'AA')
D = seqpdist(Seqs, ...'Scale', ScaleValue, ...)
D = seqpdist(Seqs, ...'GapOpen', GapOpenValue, ...)%空位罚分,默认为8
D = seqpdist(Seqs, ...'ExtendGap', ExtendGapValue, ...)
3、对序列比对的距离进行构树,
>>tree = seqlinkage(distances,'UPGMA',seqs)
%Tree = seqlinkage(Dist, Method, Names)
方法如下:
Names来描述分支点。
%'single'
Nearest distance (single linkage method)
%'complete'
Furthest distance (complete linkage method)
%'average' (default)
Unweighted Pair Group Method Average (UPGMA, group average).
%'weighted'
Weighted Pair Group Method Average (WPGMA)
%'centroid'
Unweighted Pair Group Method Centroid (UPGMC)
%'median'
Weighted Pair Group Method Centroid (WPGMC)
4、画出系统发育树
>>h = plot(tree,'orient','bottom');
>>ylabel('Evolutionary distance')
>>set(h.terminalNodeLabels,'Rotation',-45)
5、下面加上7条序列
>>data2 = {'Puti_Orangutan' 'AF451972';
'Jari_Orangutan' 'AF451964';
'Western_Lowland_Gorilla' 'AY079510';
'Eastern_Lowland_Gorilla' 'AF050738';
'Chimp_Schweinfurthii' 'AF176722';
'Chimp_Vellerosus' 'AF315498';
'Chimp_Verus' 'AF176731';
};
6、获取序列
>>for ind = 1:7
seqs(ind+5).Header = data2{ind,1};
seqs(ind+5).Sequence = getgenbank(data2{ind,2},...
'sequenceonly', true);
end
7、同上操作
>>distances = seqpdist(seqs,'Method','Jukes-Cantor','Alpha','DNA');
>>tree = seqlinkage(distances,'UPGMA',seqs);
8、同上画图
>>h = plot(tree,'orient','bottom');
>>ylabel('Evolutionary distance')
>>set(h.terminalNodeLabels,'Rotation',-45)
9、罗列系统发育树的名称
>>names = get(tree,'LeafNames')
10、选择输出
>>[h_all,h_leaves] = select(tree,'reference',3,...
'criteria','distance',...
'threshold',0.6);
11、去除无用的枝杈
>>leaves_to_prune = ~h_leaves;
pruned_tree = prune(tree,leaves_to_prune)
h = plot(pruned_tree,'orient','bottom');
ylabel('Evolutionary distance')
set(h.terminalNodeLabels,'Rotation',-30)
12、显示最终结果
phytreetool(pruned_tree)