在 Matlab 中,我想做一个 seqlogo plot氨基酸序列谱。但是,我希望所有的列都具有相同的高度,而不是通过熵来缩放绘图列的高度。
我正在修改 this question 的答案中的代码,但我想知道是否有 seqlogo 的参数或我错过的其他一些函数可以使列高统一。
或者,是否有统计转换可以应用于序列配置文件以破解所需的输出? (列高一致,每个字母的高度成线性比例
它在 seqprofile 中的概率)
最佳答案
解决这个问题的最简单方法可能是直接修改 Bioinformatics Toolbox function SEQLOGO 的代码。 (如果可能的话)。在 R2010b 中,您可以执行以下操作:
edit seqlogo
该函数的代码将显示在编辑器中。接下来,找到以下几行(第 267-284 行)并将它们注释掉或完全删除它们:
S_before = log2(nSymbols);
freqM(freqM == 0) = 1; % log2(1) = 0
% The uncertainty after the input at each position
S_after = -sum(log2(freqM).*freqM, 1);
if corrError
% The number of sequences correction factor
e_corr = (nSymbols -1)/(2* log(2) * numSeq);
R = S_before - (S_after + e_corr);
else
R = S_before - S_after;
end
nPos = (endPos - startPos) + 1;
for i =1:nPos
wtM(:, i) = wtM(:, i) * R(i);
end
然后将此行放在他们的位置:
wtM = bsxfun(@times,wtM,log2(nSymbols)./sum(wtM));
您可能希望以新名称保存文件,例如
seqlogo_norm.m
,因此您仍然可以使用未经修改的原始 SEQLOGO功能。现在您可以创建所有列归一化为相同高度的序列剖面图。例如:S = {'LSGGQRQRVAIARALAL',... %# Sample amino acid sequence
'LSGGEKQRVAIARALMN',...
'LSGGQIQRVLLARALAA',...
'LSGGERRRLEIACVLAL',...
'FSGGEKKKNELWQMLAL',...
'LSGGERRRLEIACVLAL'};
seqlogo_norm(S,'alphabet','aa'); %# Use the modified SEQLOGO function
旧答案:
我不确定如何转换序列配置文件信息以从 Bioinformatics Toolbox function SEQLOGO 获得所需的输出,但我可以向您展示如何修改替代方案
seqlogo_new.m
我为 my answer to the related question 写的你链接到。如果您更改初始化 bitValues
的行由此:bitValues = W{2};
对此:
bitValues = bsxfun(@rdivide,W{2},sum(W{2}));
然后你应该将每一列缩放到 1 的高度。例如:
S = {'ATTATAGCAAACTA',... %# Sample sequence
'AACATGCCAAAGTA',...
'ATCATGCAAAAGGA'};
seqlogo_new(S); %# After applying the above modification
关于Matlab:具有统一绘图列高度的 seqlogo,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/5133158/