系统发生分析 -开始太傻了！！

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)