MATLAB Answers

0

Generate neighbors of a string

Asked by Paolo Binetti on 9 Apr 2017
Latest activity Commented on by Paolo Binetti on 3 Jun 2017
I need to generate a cell array (or char array) containing all d-neighbors of a string of lenght k. This means all the k-long strings which differ from the input at most d characters.
My specific need concerns an alphabet of four letters, A, C, G, T. Here is an example: with string='CTA' and d=2 as inputs, the output array should have 37 elements, like 'AGA', 'CAA', 'GGA', ... and 'CTA'.
My code, attached, is the fastest of several I have tried. The main function "neighbors" calls auxiliary function "immediate_neighbors". The expected inputs of these function are provided as comments in the script.
I am convinced there are ways to achieve this much faster, but could not find any.

  3 Comments

Jan
on 9 Apr 2017
It would be much easier to understand the code if there is any documentation about the inputs. I'm not aware of a "j++" operation in Matlab, so are you sure that you've posted running code?
Some example input and output data would be useful also.
"I'm not aware of a "j++" operation in Matlab"
There isn't in MATLAB, but there is in Octave.
Thank you for the feedback. I have documented the inputs in the function script, corrected the code, and tested on Matlab.

Sign in to comment.

1 Answer

Answer by Guillaume
on 9 Apr 2017
Edited by Guillaume
on 9 Apr 2017
 Accepted Answer

For short strings and small d, your code is slightly faster, but for longer strings with larger d (tested with 18 chars and d = 3) the following is much faster:
function patterns = neighbours(pattern, distance, charset)
if nargin < 3
charset = 'ACGT';
end
%build replacement list for each character:
charsetreps = arrayfun(@(cidx) charset([1:cidx-1, cidx+1:end]), 1:numel(charset), 'UniformOutput', false);
charsetreps = vertcat(charsetreps{:});
%fill patterns for distances from 0 to distance
patterns = {pattern}; %0 distance
for d = 1:distance
%get all combinations of d indices of pattern to replace at once
charidx = nchoosek(1:numel(pattern), d);
%storage for all replacement at d:
allreps = cell(size(charidx, 1), (numel(charset)-1)^d);
%get cartesion product of replacement column in charsetreps
repcols = cell(1, d);
[repcols{:}] = ndgrid(1:numel(charset)-1);
repcols = reshape(cat(d+1, repcols{:}), [], d);
repcols = numel(charset) * (repcols-1); %conversion to linear indexing of the column
%iterate over the charidx combinations
for ci = 1:size(charidx, 1)
[~, reprows] = ismember(pattern(charidx(ci, :)), charset);
%iterate over the replacement characters
for rep = 1:size(repcols, 1)
replacement = pattern;
replacement(charidx(ci, :)) = charsetreps(reprows + repcols(rep, :)); %reprow + repcols creates linear indices
allreps{ci, rep} = replacement;
end
end
patterns = [patterns; allreps(:)]; %#ok<AGROW>
end
end
It is also a lot more generic as it works with any character set, not just 'ACGT'.
However, note that your own code could be optimised by taking the unique and assignment to f out of the loop.

  2 Comments

Guillaume, thank you for your help. Your code is indeed faster: with d = 3, I start seeing a difference when patterns are longer than 6 chars. Still, it takes quite some time for a task that seemed so simple at first sight.
As for my code, I have tried taking "unique" out of the loop, but it's actually much slower. Perhaps I did not understand your suggestion.
Guillaume, after spending a while thinking of better ways to do it, I decided to take a closer look at your code. I was able to vectorize most of the inner loop, so now the code is significantly faster even with shorter strings:
function patterns = neighbors_4(pattern, distance, charset)
if nargin < 3
charset = 'ACGT';
end
%build replacement list for each character:
charsetreps = arrayfun(@(cidx) charset([1:cidx-1, cidx+1:end]), 1:numel(charset), 'UniformOutput', false);
charsetreps = vertcat(charsetreps{:});
%fill patterns for distances from 0 to distance
patterns = {pattern}; %0 distance
for d = 1:distance
%get all combinations of d indices of pattern to replace at once
charidx = nchoosek(1:numel(pattern), d);
%storage for all replacement at d:
allreps = cell(size(charidx, 1), (numel(charset)-1)^d);
%get cartesion product of replacement column in charsetreps
repcols = cell(1, d);
[repcols{:}] = ndgrid(1:numel(charset)-1);
repcols = reshape(cat(d+1, repcols{:}), [], d);
repcols = numel(charset) * (repcols-1); %conversion to linear indexing of the column
rep_mat = pattern(ones(size(repcols, 1), 1), :);
if d == 1
[~, reprows] = ismember(pattern(charidx)', charset);
else
[~, reprows] = ismember(pattern(charidx), charset);
end
%iterate over the charidx combinations
for ci = 1:size(charidx, 1)
ccc = charidx(ci, :); %%%
replacement_mat = rep_mat;
replacement_mat(:, ccc) = charsetreps(reprows(ci, :) + repcols);
allreps(ci, :) = cellstr(replacement_mat)';
end
patterns = [patterns; allreps(:)]; %#ok<AGROW>
end
patterns = patterns';
end
In addition, I noticed that the computation of charsetreps is independent of the specific input sequence, and it is relatively time-consuming. Therefore, if one needs to call the function multiple times (for example to compute the neighbours of hundreds of different sequences), it is better to compute charsetreps just once outside of the function and provide it as an input.
So, ultimately, your function with these two changes, helped increade the speed of the main code using this function by about 10x.
Thank you again and congrats for providing the function so quickly!

Sign in to comment.