read_seqs.m

greedy em 混和模型训练算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% READ THE TEXT FILE OF SEQUENCES AND CONSRUCT %%%%%%%%%%%%%
%%%%%% THE SUBSTRINGS OF LENGTH W                   %%%%%%%%%%%%%
%%%%%% FOR DISCOVERING MOTIFS OF LENGTH W           %%%%%%%%%%%%%
%%%%%%                                                      %%%%%
%%%%%% Kostas Blekas, 19 May 2001                           %%%%%
%%%%%% please contact at [e-mail] in case of problems %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [new_seq, n, thesh, background_gen, seq, seqname, fromseq, frompos, lenmin, lenmax, alphabet] = read_seqs(W,fdata);

% alphabet ( for DNA A = {'A','C','G','T'} );
%A = {'a','c','d','e','f','g','h','i','k','l','m','n','p','q','r','s','t','v','w','y'};
A = {'A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y'};
alphabet = char(A)';
L = size(alphabet,2);

% read the text file of sequences

fid = fopen(fdata,'r');

k = 1;
while ~feof(fid)
    s1 = fscanf(fid,'%s*[^\n]');
    if s1(1) ~= '>'
        seq1(k) = cellstr(s1);
        lengthseq(k) = length(s1);
        k = k + 1;
    else
        seqname(k) = cellstr(s1);
    end
end

seq = char(seq1);

N = size(seq,1);

k1=1;
k2 = 0;
for i=1:N
    for j=1:lengthseq(i)-W+1
         new_seq(k1,:) = seq(i,j:j+W-1); 
         lenmin(k1) = k2;
         lenmax(k1) = k2 + lengthseq(i)-W+1;
         fromseq(k1) = i;
         frompos(k1) = j;
         k1 = k1+1;
     end
     k2 = k2 + lengthseq(i) - W + 1;    
 end

n=size((new_seq),1);
thesh = zeros(n,W);

for i=1:n
    for l=1:W
        thesh(i,l) = find(alphabet==new_seq(i,l));
    end
end

ll = zeros(1,L);
for i=1:L
    ll(i) = size(find(new_seq==alphabet(i)),1); 
end
background_gen = ll ./ sum(ll);

fprintf('Read N=%d sequences and create a training set of n=%d substrings of length W=%d\n',N,n,W);