Generate neighbors of a string

2 ビュー (過去 30 日間)
Paolo Binetti
Paolo Binetti 2017 年 4 月 9 日
コメント済み: Paolo Binetti 2017 年 6 月 3 日
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 件のコメント
Stephen23
Stephen23 2017 年 4 月 9 日
"I'm not aware of a "j++" operation in Matlab"
There isn't in MATLAB, but there is in Octave.
Paolo Binetti
Paolo Binetti 2017 年 4 月 9 日
編集済み: Paolo Binetti 2017 年 4 月 9 日
Thank you for the feedback. I have documented the inputs in the function script, corrected the code, and tested on Matlab.

サインインしてコメントする。

採用された回答

Guillaume
Guillaume 2017 年 4 月 9 日
編集済み: Guillaume 2017 年 4 月 9 日
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 件のコメント
Paolo Binetti
Paolo Binetti 2017 年 4 月 16 日
編集済み: Paolo Binetti 2017 年 4 月 16 日
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.
Paolo Binetti
Paolo Binetti 2017 年 6 月 3 日
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!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeGenomics and Next Generation Sequencing についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by