MATLAB Базовая реализация цепи Маркова

Question

Я пишу код, имитирующий очень простую цепь Маркова для генерации 10000 6-нуклеотидных последовательностей из любой из двух матриц перехода (т. Е. Если предыдущий нуклеотид был А, то используйте этот набор вероятностей для генерации следующего нуклеотида и т. Д.). У меня также есть запущенные продукты, которые я получаю, выбирая соответствующие вероятности из указанных матриц, но они также не распространяются правильно.

Я знаю, что MATLAB может не обрабатывать строки / символы как другие языки, поэтому я не совсем уверен, что происходит с моим кодом. По сути, я использую выборку, но, похоже, она не распределяется правильно. Мой код ниже. Это по общему признанию не элегантно, но это должно работать... но это не так. Есть ли более простой способ сделать это, учитывая эти матрицы перехода? (Очевидно, что код работает только для первой матрицы (M1).)

clear all;
close all;

%length of sequences
n = 6;

%nucleotide map

N = 'ACGT';

%initial probabilities
%     A    C    G    T
M0 = [0.25 0.25 0.25 0.25];

%transition probabilities (if row, then col)
%inside CpG
%     A       C       G       T
M1 = [0.18    0.274   0.426   0.12    %A
      0.171   0.367   0.274   0.188   %C
      0.161   0.339   0.375   0.125   %G
      0.079   0.355   0.384   0.182]; %T

%outside CpG
%     A       C       G       T
M2 = [0.3     0.205   0.285   0.21    %A
      0.322   0.298   0.078   0.302   %C
      0.248   0.246   0.298   0.208   %G
      0.177   0.239   0.292   0.292]; %T

sampletype = input('Matrix Type (M1 or M2): ', 's');

sampseq = zeros(n,1);
PM1 = 1;
PM2 = 1;
count = 0;

if strcmp(sampletype,'M1')
    for i = 1:10000
        for position = 1:n
            if position == 1
                firstnuc = randsample(N,1,true,M0);
                sampseq(1) = firstnuc;

                PM1 = PM1*0.25;
                PM2 = PM2*0.25;

            %check for previous nucleotide
            elseif position ~= 1
                if strcmp(sampseq(position-1),'A')
                    tempnuc = randsample(N,1,true,M1(1,:));
                    sampseq(position) = tempnuc;

                    PM1 = PM1*M1(1,findstr(N,tempnuc));
                    PM2 = PM2*M2(1,findstr(N,tempnuc));
                elseif strcmp(sampseq(position-1),'C')
                    tempnuc = randsample(N,1,true,M1(2,:));
                    sampseq(position) = tempnuc;

                    PM1 = PM1*M1(2,findstr(N,tempnuc));
                    PM2 = PM2*M2(2,findstr(N,tempnuc));
                elseif strcmp(sampseq(position-1),'G')
                    tempnuc = randsample(N,1,true,M1(3,:));
                    sampseq(position) = tempnuc;

                    PM1 = PM1*M1(3,findstr(N,tempnuc));
                    PM2 = PM2*M2(3,findstr(N,tempnuc));
                elseif strcmp(sampseq(position-1),'T')
                    tempnuc = randsample(N,1,true,M1(4,:));
                    sampseq(position) = tempnuc;

                    PM1 = PM1*M1(4,findstr(N,tempnuc));
                    PM2 = PM2*M2(4,findstr(N,tempnuc));
                end
            end
        end

        if PM2 > PM1
            count = count + 1;
        end
        PM1
        PM2
        PM1 = 1;
        PM2 = 1;
        sampseq
        sampseq = zeros(n,1);
    end
end

count

Answer 1

Ах, я исправил это. Чтобы обойти битовые обозначения, созданные строкой 'ACGT', я просто создал заполнители:

%    A C G T
N = [1 2 3 4]

Исправление индексации было тривиальным после этого.

Мне все еще интересно выяснить, как генерировать "реальные" последовательности, хотя, как я изначально пытался это сделать, вместо просто строк комбинаций 1, 2, 3 и 4.