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

news2025/1/18 8:43:29

这次开发了一个生物信息学比较常用的序列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_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/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

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/155613.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

再见2022，你好2023:八年程序媛老兵的践行、思考与展望

再见2022，你好2023:八年程序媛老兵的践行、思考与展望

再见2022&#xff0c;你好2023写在前面的话1.2022速记1.1 产假前&#xff0c;ParaView三维自动化项目1.2 产假后&#xff0c;EDA仿真项目1.3 从EDA行业谈谈2022年的经济寒冬2. 2023年的新年flag2.1 flag one:挑战高薪2.2 flag two:读更多的书&#xff0c;读更多专业书2.2.1 读过…
阅读更多...
动态内存管理题目讲解

动态内存管理题目讲解

前言&#xff1a; 上一期我们讲述了有关动态内存管理的知识点&#xff0c;这期我们通过几个经典的笔试题来进行深入的了解以及对知识点的巩固 目录第一题第二题第三题第四套第一题 试题如下&#xff1a; void GetMemory(char* p) {p (char*)malloc(100); } void Test(void) …
阅读更多...
1.H3CNE-计算机网络概述

1.H3CNE-计算机网络概述

计算机网络概述计算机网络定义一组自治计算机互联的集合计算机网络基本功能资源共享综合信息服务分布式处理与负载均衡计算机网络的类型局域网LAN&#xff08;Local Area Network) 由用户自行建设&#xff0c;使用私有地址组建的网络城域网MAN(Metropolitan Area Network)由运营…
阅读更多...
为什么职场第一开发语言会是SQL？看完这些你就瞬间明白了

为什么职场第一开发语言会是SQL？看完这些你就瞬间明白了

看到一个有趣的比喻&#xff0c;用来说明SQL与Excel的差别是什么。 如果把SQL比作火车&#xff0c;把Excel更比作卡车。 卡车灵活自由&#xff0c;高速或乡村小道想去哪就去哪&#xff0c;但即便每天不停歇卡车的运载量也不大&#xff0c;而且容易出交通事故。 火车运载量大…
阅读更多...
windows下软件安装：miniconda下安装R4.1.3并将其添加到 Jupyter notebook 中

windows下软件安装：miniconda下安装R4.1.3并将其添加到 Jupyter notebook 中

0. 说明&#xff1a; 本来是想在windows中用conda安装R&#xff0c;然后再下载安装RStudio并对其进行配置使之可以用conda环境中的R&#xff0c;但是经过尝试以及网络上查找相关文档发现&#xff0c;原版RStudio不支持使用conda环境中的R&#xff08;可能Anaconda中自带的RStu…
阅读更多...
Visual studio C++桌面应用程序添加外部文件引用

Visual studio C++桌面应用程序添加外部文件引用

C桌面应用程序添加外部文件引用 前言 之前对C的开发接触很少&#xff0c;本章节记录一下Visual studio开发C桌面应用程序是如何引入外部文件 ★提高阅读体验★ &#x1f449; ♠一级标题 &#x1f448; &#x1f449; ♥二级标题 &#x1f448; &#x1f449; ♥ 三级标…
阅读更多...
Apache Spark 机器学习 基本统计 1

Apache Spark 机器学习 基本统计 1

1 基本概念 相关性&#xff0c;是指两个变量或者两个系列变量的关联程度&#xff0c;也就是&#xff0c;其中一方变量的变化会影响另外一方变量的变化。 相关性分为三种关系&#xff0c;正相关、负相关以及不相关。 正相关&#xff0c;从单调递增的角度看&#xff0c;其中一…
阅读更多...
Netty基础入门——文件编程、网络编程【2】

Netty基础入门——文件编程、网络编程【2】

Netty基础入门——文件编程、网络编程【2】 基础入门【1】 1 文件编程 1.1 channel 两个channel传输数据 transferTo方法一次性最多传输2G大小的文件&#xff0c;如果超出会丢弃 public static void main(String[] args) {try (FileChannel from new FileInputStream(&quo…
阅读更多...
APM系统是什么？有什么用处？

APM系统是什么？有什么用处？

自SpringCloud问世以来&#xff0c;微服务以席卷之势风靡全球&#xff0c;企业架构都在从传统SOA向微服务转型。然而微服务这把双刃剑在带来各种优势的同时&#xff0c;也给运维、性能监控、错误的排查带来的极大的困难。在大型项目中&#xff0c;服务架构会包含数十乃至上百个…
阅读更多...
分布式助力光伏太阳能规模化发展解决方案

分布式助力光伏太阳能规模化发展解决方案

行业背景 光伏太阳能作为一种清洁环保的能源&#xff0c;得到各种开发利用&#xff0c;光伏太阳能电池板是其中的重点研究对象&#xff0c;其质量是影响太阳能电池发电效率的主要因素,所以对电池板表面质量的检测是生产中一个重要环节。随着工业的发展&#xff0c;太阳能电池板…
阅读更多...
198：vue+openlayers 解决drawend后不能获取当前feature的方法

198：vue+openlayers 解决drawend后不能获取当前feature的方法

第198个 点击查看专栏目录 本示例的目的是介绍如何在vue+openlayers项目中绘制矩形,drawend触发事件，要获取到当前绘制的feature的信息。drawend触发的时刻，add feature to the source or collection 这个变化还没有发生，所以用source.getFeatures()是获取不到最新数据的。可…
阅读更多...
OpenStack 认证Api

OpenStack 认证Api

在调用OpenStack的Api或者其它组建的Api时都需要进行 OpenStack 认证&#xff0c;在这里记录一下如何调用OpenStack 认证接口或者token 和给其它接口增加token的方式一. 调用OpenStack auth接口接口地址&#xff1a;http://ip:5000/v3/auth/tokens参数&#xff1a;{"auth&…
阅读更多...
特色风情小镇行业发展动态及市场需求前景分析

特色风情小镇行业发展动态及市场需求前景分析

2023-2029年中国特色风情小镇行业发展动态及市场需求前景报告报告编号&#xff1a;1691653免费目录下载&#xff1a;http://www.cninfo360.com/yjbg/qthy/qt/20230110/1691653.html本报告著作权归博研咨询所有&#xff0c;未经书面许可&#xff0c;任何组织和个人不得以任何形式…
阅读更多...
结构体内存对齐与结构体位段：学习笔记8

结构体内存对齐与结构体位段：学习笔记8

目录 一.结构体基础知识 1. 结构体的特殊声明 2. 结构的自引用 3.结构体变量的定义和初始化 二.结构体内存对齐 1.关键概念&#xff1a; 2.计算示例 3.嵌套结构体的内存计算 4.结构体内存对齐的意义 5.定义结构体时的注意事项 6.修改默认对齐数 附&#xff1a;关…
阅读更多...
【PWA学习】5. 使用 Notification API 来进行消息提醒

【PWA学习】5. 使用 Notification API 来进行消息提醒

引言 在上一节, 介绍了如何使用 Push API 进行服务端消息推送。提到 Push 就不得不说与其联系紧密的另一个 API——Notification API。它让我们可以在“网站外”显示消息提示&#xff1a; 消息推送示例即使当你切换到其他 Tab&#xff0c;也可以通过提醒交互来快速让用户回到你…
阅读更多...
webviz安装，docker安装可正常使用与Foxglove Studio

webviz安装，docker安装可正常使用与Foxglove Studio

Foxglove Studio Foxglove Studio与webviz使用起来非常类似 去可以直接使用web也可以下载安装包 Foxglove Studio不提供源码 安装包下载地
阅读更多...
linux cgroup、kubernetes limit

linux cgroup、kubernetes limit

linux cgroup、kubernetes limit 1.cgroups 简介 cgroups&#xff0c;其名称源自控制组群&#xff08;control groups&#xff09;的缩写&#xff0c;是内核的一个特性&#xff0c;用于限制、记录和隔离一组进程的资源使用&#xff08;CPU、内存、磁盘 I/O、网络等&#xff0…
阅读更多...
JSP——分页查询

JSP——分页查询

✅作者简介&#xff1a;热爱国学的Java后端开发者&#xff0c;修心和技术同步精进。 &#x1f34e;个人主页&#xff1a;Java Fans的博客 &#x1f34a;个人信条&#xff1a;不迁怒&#xff0c;不贰过。小知识&#xff0c;大智慧。 &#x1f49e;当前专栏&#xff1a;JAVA开发者…
阅读更多...
Homekit智能家居DIY产品一智能面板开关

Homekit智能家居DIY产品一智能面板开关

触摸开关&#xff0c;即通过触摸方式控制的墙壁开关&#xff0c;其感官场景如同我们的触屏手机&#xff0c;只需手指轻轻一点即可达到控制电器的目的&#xff0c;随着人们生活品质的提高&#xff0c;触摸开关将逐渐将换代传统机械按键开关。 触摸开关控制原理 触摸开关我们把…
阅读更多...
【广度优先搜索遍历 BFS】单词接龙

【广度优先搜索遍历 BFS】单词接龙

一、题目描述 字典 wordList 中从单词 beginWord 和 endWord 的 转换序列 是一个按下述规格形成的序列 beginWord -> s1 -> s2 -> ... -> sk&#xff1a; - 每一对相邻的单词只差一个字母。 - 对于 1 < i < k 时&#xff0c;每个 si 都在 wordList 中。注意…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023