这次开发了一个生物信息学比较常用的序列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文件中读取数据，若数据是这样一行名称一行序列的格式：

>seq01ACCCTGTAAGTTTT>seq02TCAGTGTAAGTATC>seq03CATTCGTAAGTACC>seq04CGCTGGTAAGGACT>seq05ACCGGGTGAGCGCG

则可通过如下代码读取：

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

序列txt

若txt文件中只有序列：

GNYLGLTVETISRLLASYLGLRLETVCRSVSEMTGTTLHTVSRLLAEMTGTTLHTVSRILAEMTGTTLHTVSRILAEMTGTTVETTIRVMASRVGLTVQTVSTIVAARLGLTPETFSRVLADYLGTTPETVSRTLADMLGSKRETVSRQLANYIGTSPETISRKIATFIGTTPETISRKFSAFIGTTPETISRKLADVLGLSVVHMNRVIADALGLTPIHINRMLAEAIGSTRVTVTRLLGNYLGLTVETISRLLGNYLGLTVETISRLLGNYLGLTIETISRLLGNYLGLTIETISRLLGNYLGLTVETISRLLGNYLGLTVETISRLL

则可通过如下代码读取：

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_n\right) 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/PRcoe.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};    endend% 获取版本信息tver=version('-release');verMatlab=str2double(tver(1:4))+(abs(tver(5))-abs('a'))/2;if verMatlab<2017    hold onelse    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};    endend% 统计基因出现次数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)];    endendcoe.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')    endend% 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