MATLAB | 如何使用MATLAB绘制序列logo图

admin • • 5G

这次开发了一个生物信息学比较常用的序列logo图绘制MATLAB代码包，绘制效果如下：

数据来自基迪奥生物项目编号为seqlogojrois9l2jit的示例数据。同时本工具函数参考以下文献：

  • Tareen A, Kinney J B. Logomaker: beautiful sequence logos in Python[J]. Bioinformatics, 2020, 36(7): 2272-2274.
  • Crooks G E, Hon G, Chandonia J M, et al. WebLogo: a sequence logo generator[J]. Genome research, 2004, 14(6): 1188-1190.

教程部分

1 数据格式

char矩阵

数据使用char类型矩阵，每一行代表一条序列：

Data =['GNYLGLTVETISRLL';
    'GNYLGLTVETISRLL';
    'GNYLGLTVETISRLL';
    'GNYLGLTIETISRLL']

名称与序列txt

当然我们可以从txt文件中读取数据，若数据是这样一行名称一行序列的格式：

>seq01
ACCCTGTAAGTTTT
>seq02
TCAGTGTAAGTATC
>seq03
CATTCGTAAGTACC
>seq04
CGCTGGTAAGGACT
>seq05
ACCGGGTGAGCGCG

则可通过如下代码读取：

Data=readcell('seqlogo_DNA.txt');
Data=Data(2:2:end,1);
Data=reshape([Data{:}],[],length(Data))';

序列txt

若txt文件中只有序列：

GNYLGLTVETISRLL
ASYLGLRLETVCRSV
SEMTGTTLHTVSRLL
AEMTGTTLHTVSRIL
AEMTGTTLHTVSRIL
AEMTGTTVETTIRVM
ASRVGLTVQTVSTIV
AARLGLTPETFSRVL
ADYLGTTPETVSRTL
ADMLGSKRETVSRQL
ANYIGTSPETISRKI
ATFIGTTPETISRKF
SAFIGTTPETISRKL
ADVLGLSVVHMNRVI
ADALGLTPIHINRML
AEAIGSTRVTVTRLL
GNYLGLTVETISRLL
GNYLGLTVETISRLL
GNYLGLTIETISRLL
GNYLGLTIETISRLL
GNYLGLTVETISRLL
GNYLGLTVETISRLL

则可通过如下代码读取：

Data=readcell('seqlogo_protein_3.txt');
Data=reshape([Data{:}],[],length(Data))';

2 单位显示

有Bits和Probability供用户选择，默认bits。logo图纵坐标的单位常见有两种，一种是百分比，另一种是Bits。对于Probability，很好理解，每个字母的出现频率；对于Bits，可参考下面的公式：

R

s

e

q

=

S

max 

−

S

o

b

s

=

log

⁡

2

N

−

(

−

∑

n

=

1

N

p

n

log

⁡

2

p

n

)

R_{s e q}=S_{text {max }}-S_{o b s}=log _2 N-left(-sum_{n=1}^N p_n log _2 p_nright)

Rseq​=Smax ​−Sobs​=log2​N−(−n=1∑N​pn​log2​pn​)

p

n

p_n

pn​是相应位置n上相应字符出现频率，N是不同字符的总数量（核酸为4，蛋白质为20）。因此，对于图中的y轴的最大数值就不难理解，核酸序列的最大值为

log

⁡

2

4

=

2

log_2 4 = 2

log2​4=2bits，蛋白序列为

log

⁡

2

20

≈

4.32

log_2 20≈4.32

log2​20≈4.32 bits。

可以通过设置函数的Method属性来调整显示单位，可选值为'Bits'/'Prob',其中Bits为默认值。

3 核酸序列绘制示例

将Type属性设置为NA即可绘制核酸序列logo图，当然参数的默认值就是NA所以不设置也可以，以下给出绘制核酸序列显示单位分别为Bits和Probability的绘制代码和绘图结果：

Data=readcell('seqlogo_DNA.txt');
Data=Data(2:2:end,1);
Data=reshape([Data{:}],[],length(Data))';

figure()
seqLogo(Data)

figure()
seqLogo(Data,'Method','Prob')

4 蛋白序列绘制示例

将Type属性设置为PR即可绘制蛋白序列logo图：

Data=readcell('seqlogo_protein_2.txt'); 
Data=Data(2:2:end,1);
Data=reshape([Data{:}],[],length(Data))';

figure()
seqLogo(Data,'Type','PR')

figure()
seqLogo(Data,'Type','PR','Method','Prob')

5 核酸序列配色

对于每个字母都可以单独设置颜色，自由度比较高，举个例子，将颜色设置为:
{‘C’,[205,255,101]./255,
‘A’,[104,101,255]./255,
‘TU’,[164,230,101]./255,
‘G’,[104,203,254]./255}：

Data=readcell('seqlogo_DNA.txt');
Data=Data(2:2:end,1);
Data=reshape([Data{:}],[],length(Data))';

Color={'C',[205,255,101]./255,'A',[104,101,255]./255,'TU',[164,230,101]./255,'G',[104,203,254]./255};
% Color={'C',[127,91,93]./255,'A',[187,128,110]./255,'TU',[197,173,143]./255,'G',[59,71,111]./255};

figure()
seqLogo(Data,'Color',Color)

figure()
seqLogo(Data,'Method','Prob','Color',Color)

设置为：
{‘C’,[127,91,93]./255,
‘A’,[187,128,110]./255,
‘TU’,[197,173,143]./255,
‘G’,[59,71,111]./255}

5 蛋白序列配色

对于每个字母都可以单独设置颜色，自由度比较高，下面可以设置任意多组颜色，随意划分分组：

Data=readcell('seqlogo_protein_3.txt');
Data=reshape([Data{:}],[],length(Data))';

Color={'CDEFH',[205,255,101]./255,'AIKLM',[104,101,255]./255,'TUNPQR',[164,230,101]./255,'GSVWY',[104,203,254]./255};
% Color={'CDEFH',[127,91,93]./255,'AIKLM',[187,128,110]./255,'TUNPQR',[197,173,143]./255,'GSVWY',[59,71,111]./255};

figure()
seqLogo(Data,'Type','PR','Color',Color)

figure()
seqLogo(Data,'Type','PR','Method','Prob','Color',Color)

.

6 子图

当然subplot函数啥的也可以用：

Data=readcell('seqlogo_DNA.txt');
Data=Data(2:2:end,1);
Data=reshape([Data{:}],[],length(Data))';

figure()
subplot(2,1,1)
seqLogo(Data)

subplot(2,1,2)
seqLogo(Data,'Method','Prob')

工具函数完整代码

仅代码无法运行，需要文件夹内有logoData.mat文件，此处仅做展示，完整代码+mat文件+示例数据可以从文末提供的MATHWORKS的fileexchange链接地址下载，或者从文末网盘链接下载。

function seqLogo(varargin)
% @author : slandarer
% gzh  : slandarer随笔
% Zhaoxu Liu / slandarer (2023). sequence logos (序列logo图) 
% (https://www.mathworks.com/matlabcentral/fileexchange/123060-sequence-logos-logo), 
% MATLAB Central File Exchange. 检索来源 2023/1/10.
% =========================================================================
% Color  | 配色              | {'A',[]./255,'C',[]./255,... ...}
% Method | 比例计算方法       | Bits/Prob
% Type   | 种类(核酸/蛋白质)  | NA/PR
coe.arginList={'Color','Method','Type'};
% 数据预定义
logoData=load('logoData.mat');
coe.Color={'CDEFH',[37,92,153]./255,'AIKLM',[16,150,72]./255,'TUNPQR',[214,40,57]./255,'GSVWY',[247,179,43]./255};
coe.Method='Bits';
coe.Type='NA';
% 坐标区域获取
if isa(varargin{1},'matlab.graphics.axis.Axes')
    coe.ax=varargin{1};varargin(1)=[];
else
    coe.ax=gca;
end
% 获取其他数据
coe.Data=varargin{1};varargin(1)=[];
for i=1:2:(length(varargin)-1)
    tid=ismember(coe.arginList,varargin{i});
    if any(tid)
        coe.(coe.arginList{tid})=varargin{i+1};
    end
end
% 获取版本信息
tver=version('-release');
verMatlab=str2double(tver(1:4))+(abs(tver(5))-abs('a'))/2;
if verMatlab<2017
    hold on
else
    hold(coe.ax,'on')
end
% 颜色计算
coe.CData=zeros(length(logoData.logoName),3);
for i=1:2:length(coe.Color)
    tLogo=coe.Color{i};
    for j=1:length(tLogo)
        tPos=find(logoData.logoName==tLogo(j));
        coe.CData(tPos,:)=coe.Color{i+1};
    end
end
% 统计基因出现次数
coe.Count=zeros(length(logoData.logoName),size(coe.Data,2));
for i=1:length(logoData.logoName)
    coe.Count(i,:)=sum(coe.Data==logoData.logoName(i),1);
end
% 开始绘图
if strcmpi(coe.Method,'Prob')
    coe.ax.YLim=[0,1];
    coe.ax.DataAspectRatio=[1 .2 1];
    coe.ax.YLabel.String='Probability';
else
    coe.ax.YLabel.String='Bits';
    if strcmpi(coe.Type,'NA')
        coe.ax.DataAspectRatio=[1 .4 1];
        coe.ax.YLim=[0,log(4)/log(2)];
    else
        coe.ax.DataAspectRatio=[1 .8 1];
        coe.ax.YLim=[0,log(20)/log(2)];
    end
end
coe.ax.XLim=[.5,size(coe.Data,2)+.5];
coe.ax.XTick=1:size(coe.Data,2);
coe.ax.LineWidth=1.2;
coe.ax.TickDir='out';
coe.ax.TickLength=[0.0020 0.0250];
coe.ax.FontSize=14;
coe.ax.YLabel.FontSize=16;
%coe.ax.LooseInset=[0,0,0,0];
fig=gcf;
fig.Units='normalized';
fig.Position=[0,0,1,1];
for i=1:size(coe.Count,2)
    tPos=find(coe.Count(:,i)>0);
    tCount=coe.Count(tPos,i)';
    tRatio=tCount./sum(tCount);
    if strcmpi(coe.Method,'Prob')
        maxH=1;
    else
        if strcmpi(coe.Type,'NA')
            maxH=log(4)/log(2);
        else
            maxH=log(20)/log(2);
        end
        maxH=maxH+sum(log(tRatio)./log(2).*tRatio);
    end
    tLen=tRatio.*maxH;
    [sortRatio,ind]=sort(tLen);
    cumsumSortRatio=[0,cumsum(sortRatio)];
    for j=1:length(sortRatio)
        tPic=logoData.logoPic.(logoData.logoName(tPos(ind(j))));
        tAlpha=double(255-tPic(:,:,1)./3-tPic(:,:,2)./3-tPic(:,:,3)./3)./255;
        tPic(:,:,1)=(255-tPic(:,:,1)).*coe.CData(tPos(ind(j)),1);
        tPic(:,:,2)=(255-tPic(:,:,2)).*coe.CData(tPos(ind(j)),2);
        tPic(:,:,3)=(255-tPic(:,:,3)).*coe.CData(tPos(ind(j)),3);
        image([i-.5,i+.5],[cumsumSortRatio(j),cumsumSortRatio(j+1)],tPic,'AlphaData',tAlpha,'AlphaDataMapping','scaled')
    end
end
% Zhaoxu Liu / slandarer (2023). sequence logos (序列logo图) 
% (https://www.mathworks.com/matlabcentral/fileexchange/123060-sequence-logos-logo), 
% MATLAB Central File Exchange. 检索来源 2023/1/10.
end

完

未经允许本代码请勿作商业用途，引用的话可以引用我file exchange上的链接，可使用如下格式：

Zhaoxu Liu / slandarer (2023). sequence logos (序列logo图) (https://www.mathworks.com/matlabcentral/fileexchange/123060-sequence-logos-logo), MATLAB Central File Exchange. 检索来源 2023/1/10.

若转载请保留以上file exchange链接及本文链接！！！！!

链接：
https://pan.baidu.com/s/1RMafcUm0NVyz6I0R93-DUw?pwd=slan
提取码：slan

本图文内容来源于网友网络收集整理提供，作为学习参考使用，版权属于原作者。
THE END
分享
二维码
< <上一篇
下一篇>>