Matlab: seqlogo с одинаковой высотой столбца графика

Question

В Matlab я хочу сделать seqlogo график профиля аминокислотной последовательности.Но вместо того, чтобы масштабировать высоту столбцов графика энтропией, я хочу, чтобы все столбцы были одинаковой высоты.

Я нахожусь в процессе изменения кода из ответов на этот вопрос , но мне интересно, есть ли параметр для seqlogo или какой-либо другой функции, которую я пропустил, которая сделает высоту столбцов равномерной.

В качестве альтернативы, есть ли статистическое преобразование, которое можно применить к профилю последовательностивзломать желаемый вывод?(высота столбцов одинаковая, высота каждой буквы линейно пропорциональна ее вероятности в seqprofile)

Answer 1

Вероятно, самый простой способ обойти эту проблему - это напрямую изменить код для функции Биоинформатика Toolbox 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

СТАРЫЙ ОТВЕТ:

Я не уверен, как преобразовать информацию профиля профиля, чтобы получитьжелаемый вывод из функции набора инструментов биоинформатики SEQLOGO , но я могу показать вам, как изменить альтернативу seqlogo_new.m, которую я написал для моего ответа на связанный вопрос , с которым вы связались.Если вы измените строку, которая инициализирует 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

Answer 2

На данный момент мой обходной путь - создать группу поддельных последовательностей, которые соответствуют профилю последовательности, а затем передать эти последовательности в http://weblogo.berkeley.edu/logo.cgi.Вот код для создания поддельных последовательностей:

function flatFakeSeqsFromPwm(pwm, letterOrder, nSeqsToGen, outFilename)
%translates a pwm into a bunch of fake seqs with the same probabilities
%for use with http://weblogo.berkeley.edu/

%pwm should be a 4xn or a 20xn position weight matrix. Each col must sum to 1
%letterOrder = e.g. 'ARNDCQEGHILKMFPSTWYV' for my data
%nSeqsToGen should be >= the # of pixels tall you plan to make your chart

[height windowWidth] = size(pwm);
assert(height == length(letterOrder));
assert(isequal(abs(1-sum(pwm)) < 1.0e-10, ones(1, windowWidth))); %assert all cols of pwm sum to 1.0

fd = fopen(outFilename, 'w');

for i = 0:nSeqsToGen-1
    for seqPos = 1:windowWidth
        acc = 0; %accumulator
        idx = 0;
        while i/nSeqsToGen >= acc
            idx = idx + 1;
            acc = acc + pwm(idx, seqPos);
        end
        fprintf(fd, '%s', letterOrder(idx));
    end
    fprintf(fd, '\n');
end

fclose(fd);
end