How can I find the complement of a subset of a DNA sequence using a logical index?

Question

I have a DNA sequence, its length for example is m*4n:

B = 'GATTAACTACACTTGAGGCT...';

I have also a vector of real numbers X = {xi, i = 1..m*4n}, and use mod(X,1) to keep them in the range [0,1]. For example:

X = [0.223 0.33 0.71 0.44 0.91 0.32 0.11 ....... m*4n];

I then need to transform X into a binary vector by applying the function:

f(x)={0  ,0 < X(i,j) ≤ 0.5;  1 ,0.5 < X(i,j) ≤ 1;)

The output according the previous values will be like X = [0010100 ....]. If X(i,j)==1, then B(i,j) is complemented, otherwise it is unchanged. In this context, the complement is the matching base-pair (i.e. A->T, C->G, G->C, and T->A).

This is the code I've tried so far, with no success:

%%maping X chaotic sequence from real numbers to binary sequence using threshold function
 X = v(:,3); 
 X(257)=[];
 disp (X);
 mode (X,1);
 for i=1
    for j=1:256
 if ((X(i,j)> 0) && (X(i,j)<= .5))
     X(i,j) = 0;
 elseif ((X(i,j)> .5) && (X(i,j)<= 1)) 
     X(i,j) = 1;
 end
    end
 end
 disp(X);

How can I properly perform the indexing and complement?

You'll need to be more specific, I don't even see how [image-processing] comes into this, and I don't know what you want. And even then, nobody will write your code for you, so you'll need to have a working/non-working example of some kind about which you ask a question. Otherwise your question is too broad. — Andras Deak -- Слава Україні, Dec 21 '16 at 18:31
I join Andras' comment. Please at least show a small numeric example of what you're trying to do. It's completely unclear to me what you're asking. — Dev-iL, Dec 21 '16 at 18:32
@AndrasDeak , sir, i definitely tried to write this by myself and when i failed i came here to ask the experts in this field to help me in this specific part. this is a tiny part from my whole program. thanks for help — M.A.Fathy, Dec 21 '16 at 18:58
You should add what you tried, and describe how it failed. As I said, we can help you debug your code, but we won't write it for you. If this tiny part is too hard, I'm not sure you'll have more luck with the rest. Don't you agree? — Andras Deak -- Слава Україні, Dec 21 '16 at 19:12
I'm not even sure how you're encoding nine 4-character strings into a 3x3 matrix. I can figure out your threshold function (although it would be nice to know what happens for `Z(i,j)>1`), but I'd prefer you tell us what the complement of a DNA matrix is rather than having to look it up. Sample input and desired output are essential. — beaker, Dec 21 '16 at 19:37
@M.A.Fathy: I tried to edit to improve the question. With your added code attempt, I think it deserves reopening now. If it gets reopened, I can post an answer for you. — gnovice, Dec 22 '16 at 05:49
@gnovice: thank you sir for your helping , i appreciate that! waiting for your answer . — M.A.Fathy, Dec 22 '16 at 21:17
shall i wait till it gets reopened or trying to post it again!? @gnovice — M.A.Fathy, Dec 22 '16 at 21:28
@M.A.Fathy: I'd wait. It needs one more reopen vote. If it doesn't reopen by tonight, email me and I'll send you an answer. — gnovice, Dec 22 '16 at 21:39

gnovice · Accepted Answer · 2016-12-23T04:46:47.227

Given a sample base-pair sequence stored as a character array:

B = 'GATTAACT';

And a sample vector of numeric values the same length as B:

X = [0.223 0.33 0.71 0.44 0.91 0.32 0.11 1.6];

Then there is a fairly straightforward solution...

First, your use of the mod function implies you want to use only the fractional part of each value in X. This is how you wold do that:

>> X = mod(X, 1)
X =
    0.2230    0.3300    0.7100    0.4400    0.9100    0.3200    0.1100    0.6000

Next, you should give the documentation on vectorization a read. It will teach you that a for loop can be avoided for many operations in MATLAB. In particular, applying a logical test to your vector X can be done like so:

>> index = (X > 0.5)
index =
    0   0   1   0   1   0   0   1

And index is now a logical index the same length as X with ones (i.e. true) for each value greater than 0.5. You now want to get the characters corresponding to those indices in B, change them to their complement, then place them back in B. You can do this using a little trick in MATLAB whereby a character is converted to its ASCII numeric value when used an as index:

>> compMap = '';  % Initialize to an empty string
>> compMap('ACGT') = 'TGCA'
compMap =
                                                                T G   C            A

Notice the characters 'TGCA' get placed in indices 65, 67, 71, and 84 of compMap (i.e. the ASCII values for 'ACGT'). The rest are blanks. Now, you can replace the indexed base-pairs with their complements by simply doing this:

>> B(index) = compMap(B(index))
B =
GAATTACA

Putting this all together, here's the solution:

B = '...';     % Whatever your sequence is
X = [...];     % Whatever your values are
compMap = '';
compMap('ACGT') = 'TGCA';      % Build a complement map
index = (mod(X, 1) > 0.5);     % Get your logical index
B(index) = compMap(B(index));  % Replace with complements

thanks sir for your response .. what i should change if i want to apply this on a whole matrix ( if B was a matrix not a subset of a DNA sequence) @gnovice — M.A.Fathy, Dec 23 '16 at 19:13
@M.A.Fathy: If `X` is a matrix that is the same exact size as `B`, then the code above will still work, no changes needed. — gnovice, Dec 23 '16 at 23:05
in fact, X is not the same exact size as B ... its size is 256*1 ,,, how can i make its size be 256*4(256) , is there any solution ?? i need take each 4bits of X: e.g (1001) to complement each pixel of B e.g (ACGT).. apply that for whole the matrix B ... any help plz! @gnovice — M.A.Fathy, Dec 24 '16 at 14:52

How can I find the complement of a subset of a DNA sequence using a logical index?

1 Answers1