【转】matlab学习--bioinformatics toolbox学习

转自：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)