Open main menu

Biolecture.org β

DNA Sequence Alignment

% <<<MATLAB code I created for DNA Sequence Alignment.>>>

clear; clc;
% Gene Sequence Alignment 
% 20141581_Kyunghyun Cho

fprintf('Type the sequqence of reference gene in '); cprintf('_red', 'line 10\n\n');
fprintf('Type the sequqence of your sample gene in '); cprintf('_red', 'line 13\n\n');
fprintf('Set the range of quality you are focusing in '); cprintf('_red', 'BOTH line 69, 185\n\n\n\n\n');

% <<<Type the sequence of "reference gene" here!>>>
ref = nt2int('ATGTTTGGCA'); 

% <<<Type the sequence of your "sample gene" here!>>>
sam = nt2int('ATGTC'); 


% Alignment Code - Do not edit!
[a,ref_length] = size(ref); % Basically, a=1, and the length of reference gene sequence is "ref_length"
[b,sam_length] = size(sam); % Basically, b=1, and the length of sample gene sequence is "sam_length"

% I considered the cases which the lengths of them are not identical.
% I needed some bases that mean nothing(i,e. just gap) in order to make sequences to have same length.
% Then, it is able to align two sequences that have different length.

if ref_length<sam_length
    
    length_difference = sam_length - ref_length;
    nobase = zeros(1,length_difference); % Basically, 0 means gap in sequence, so I created (1Xlength_differenece) matrix with zero.
    
    mod_ref = horzcat(ref,nobase); % Using function "horzcat" to concatenate ref and nobase marix.
    % For instance, if ref = [2,3,4] and sam = [1,2,2,3,4], then the length_difference is 2. Thus nobase = [0,0] and mod_ref = [2,3,4,0,0]. Now, size is same, then it is able to align!
    mod_sam = sam;

    % Now, then the question is how to align the sequences.
    % I used a trick using "for" function.
    
    % Using the first "for" to move the frame of "mod_sam "sequences to consider all possibility of alignment.
    % Simply, for example, suppose that our sequences are mod_ref = [2,3,4,0,0] and mod_sam = [1,2,2,3,4].
    % If just align them without consideration of possibility of frame moving, they are not matching at all. Because mod_ref(1) ~= mod_sam(1), mod_ref(2) ~= mod_sam(2) and so on.
    % However, what if the frame of mod_sam moved to like [2,3,4,1,2]. Then, three components are matching! Therefore, I considered all possible frame movement of "mod_sam".
    
    % Using the second "for" to do alignment and to calculate the quality of the alignment.
    % Simply, I created another function "aligned" which composed of all 1 or 0 values.
    % I compared two sequences whether they are matcing at certain point from begining(5') to end(3'). 1 means two sequences are matching at that point. 0 means not matcing. 
    % For example, if mod_ref = [1,2,1,4,0] and move_sam = [1,1,3,4,2], then aligned = [1,0,0,1,0].
    
    % The alignment quality I defined would be summation of all components of matrix "aligned" divide by "sam_length" multiply by 100(%).
    % I created "result" matrix that collects all possible alignment quaility during frame moving.
    % In "result", the highest value is "highest_q" which means the highest quality. And using "find" function to find where the highest_q is located in "result".

    % Using the third "for" to get most possible alignment result which has the highest quality.
    for i = 1:sam_length
        move_sam = circshift(mod_sam, i-1); 
        
        for j = 1:sam_length
            if mod_ref(j) == move_sam(j)
                aligned(j) = 1;
            else
                aligned(j) = 0;
            end
            quality = sum(aligned)/sam_length;
        end
        quality_family(i) = quality;
        
        
        
        % <<<Set the range of quality you are intersted in!>>> 
        X = find(quality_family >= 0.2);
        
        
        
        
        % Code for displaying alignment result (Aligned sequences in the range) - Do not edit!
        Number_quality_in_range = numel(X);
        
        for x = 1:Number_quality_in_range
            pos_sam_in_range(Number_quality_in_range,:) = circshift(sam, X(1,Number_quality_in_range)-1);
            
            for y = 1:sam_length
                if mod_ref(y) == pos_sam_in_range(x,y)
                pos_aligned(x,y) = 1;
                else
                pos_aligned(x,y) = 0;
                end  
            end
            
            for z = 1:sam_length
                    if pos_aligned(x,z) == 0 
                        aligned_ref.name = int2nt(mod_ref(z));
                        aligned_sam.name = int2nt(pos_sam_in_range(x,z));
                        aligned_ref.Color = [255 0 0]; 
                        aligned_sam.Color = [255 0 0];
                        pre1 = '<HTML><FONT color="';
                        post1 = '</FONT></HTML>';
                        listboxStr_ref1 = cell(numel( aligned_ref ),1);
                        listboxStr_sam1 = cell(numel( aligned_sam ),1);
                        str_ref1 = [pre1 reshape( dec2hex(aligned_ref.Color, 2)',1, 6) '">' aligned_ref.name post1];
                        listboxStr_ref1 = str_ref1;
                        Show_aligned_ref = {listboxStr_ref1};
                        str_sam1 = [pre1 reshape( dec2hex(aligned_sam.Color, 2)',1, 6) '">' aligned_sam.name post1];
                        listboxStr_sam1 = str_sam1;
                        Show_aligned_sam = {listboxStr_sam1};
                        
                    else
                        Show_aligned_ref = int2nt(mod_ref(z));
                        Show_aligned_ref = cellstr(Show_aligned_ref); 
                        Show_aligned_sam = int2nt(pos_sam_in_range(x,z));
                        Show_aligned_sam = cellstr(Show_aligned_sam);
                    end
                    Result_in_range(2*x-1,z) = Show_aligned_ref;
                    Result_in_range(2*x,z) = Show_aligned_sam; 
                end 
        end
            
        % Find highest quality - Do not edit!
        highest_q = max(quality_family); 
        [l, location_highest_q] = find(quality_family==highest_q); % Basically, l=1, so focusing on location_highest_q only.

        most_pos_sam = circshift(sam, location_highest_q-1); % Find most possible sequence of sample.
        
        for k = 1:sam_length
            if mod_ref(k) == most_pos_sam(k)
                most_aligned(k) = 1;
            else
                most_aligned(k) = 0;
            end
        end
    end
    
    % Code for displaying alignment result (Aligned sequences that have highest quality) - Do not edit!
    for m =1:sam_length
        
        if most_aligned(m) == 0 % If two sequences are different at same position, display sequences as red colour using "rgb2Hex" function.
            most_aligned_ref.name = int2nt(mod_ref(m));
            most_aligned_sam.name = int2nt(most_pos_sam(m));
            most_aligned_ref.Color = [255 0 0]; 
            most_aligned_sam.Color = [255 0 0];
            pre = '<HTML><FONT color="';
            post = '</FONT></HTML>';
            listboxStr_ref = cell(numel( most_aligned_ref ),1);
            listboxStr_sam = cell(numel( most_aligned_sam ),1);
            
            str_ref = [pre reshape( dec2hex(most_aligned_ref.Color, 2)',1, 6) '">' most_aligned_ref.name post];
            listboxStr_ref = str_ref;
            Show_most_aligned_ref = {listboxStr_ref};
            str_sam = [pre reshape( dec2hex(most_aligned_sam.Color, 2)',1, 6) '">' most_aligned_sam.name post];
            listboxStr_sam = str_sam;
            Show_most_aligned_sam = {listboxStr_sam};
            
        else
            Show_most_aligned_ref = int2nt(mod_ref(m));
            Show_most_aligned_ref = cellstr(Show_most_aligned_ref);
            Show_most_aligned_sam = int2nt(most_pos_sam(m));
            Show_most_aligned_sam = cellstr(Show_most_aligned_sam);
        end
        
    Result_highest_quality(1,m) = Show_most_aligned_ref;
    Result_highest_quality(2,m) = Show_most_aligned_sam;
    
    end    
    
    
% In the same way, but consider the case of same length or larger reference length.
else
    
    length_difference = ref_length - sam_length;
    nobase = zeros(1,length_difference);
    mod_ref = ref;
    mod_sam = horzcat(sam,nobase);
    
    for i = 1:ref_length
    move_sam = circshift(mod_sam, i-1); 
        for j = 1:ref_length
            if mod_ref(j) == move_sam(j)
                aligned(j) = 1; 
            else
                aligned(j) = 0; 
            end
        end
        quality = sum(aligned)/ref_length;
        quality_family(i) = quality; 
        
        % <<<Set the range of quality you are intersted in!>>> 
        X = find(quality_family >= 0.2);
        
        
        
        % Code for displaying alignment result (Aligned sequences in the range) - Do not edit!
        Number_quality_in_range = numel(X);

        
        for x = 1:Number_quality_in_range
            pos_sam_in_range(Number_quality_in_range,:) = circshift(mod_sam, X(1,Number_quality_in_range)-1);
            
            for y = 1:ref_length
                if mod_ref(y) == pos_sam_in_range(x,y)
                pos_aligned(x,y) = 1;
                else
                pos_aligned(x,y) = 0;
                end  
            end
            
            for z = 1:ref_length
                    if pos_aligned(x,z) == 0 
                        aligned_ref.name = int2nt(mod_ref(z));
                        aligned_sam.name = int2nt(pos_sam_in_range(x,z));
                        aligned_ref.Color = [255 0 0]; 
                        aligned_sam.Color = [255 0 0];
                        pre1 = '<HTML><FONT color="';
                        post1 = '</FONT></HTML>';
                        listboxStr_ref1 = cell(numel( aligned_ref ),1);
                        listboxStr_sam1 = cell(numel( aligned_sam ),1);
                        str_ref1 = [pre1 reshape( dec2hex(aligned_ref.Color, 2)',1, 6) '">' aligned_ref.name post1];
                        listboxStr_ref1 = str_ref1;
                        Show_aligned_ref = {listboxStr_ref1};
                        str_sam1 = [pre1 reshape( dec2hex(aligned_sam.Color, 2)',1, 6) '">' aligned_sam.name post1];
                        listboxStr_sam1 = str_sam1;
                        Show_aligned_sam = {listboxStr_sam1};
                        
                    else
                        Show_aligned_ref = int2nt(mod_ref(z));
                        Show_aligned_ref = cellstr(Show_aligned_ref); 
                        Show_aligned_sam = int2nt(pos_sam_in_range(x,z));
                        Show_aligned_sam = cellstr(Show_aligned_sam);
                    end
                    Result_in_range(2*x-1,z) = Show_aligned_ref;
                    Result_in_range(2*x,z) = Show_aligned_sam;
                end 
        end
            

        highest_q = max(quality_family);
        [l, location_highest_q] = find(quality_family==highest_q); 

        most_pos_sam = circshift(mod_sam, location_highest_q-1);
        
        for k = 1:ref_length
            if mod_ref(k) == most_pos_sam(k)
                most_aligned(k) = 1; 
            else
                most_aligned(k) = 0; 
            end
        end
    end
% Code for displaying alignment result (Aligned sequences that have highest quality) - Do not edit!
    for n =1:ref_length
        
        if most_aligned(n) == 0
            most_aligned_ref.name = int2nt(mod_ref(n));
            most_aligned_sam.name = int2nt(most_pos_sam(n));
            most_aligned_ref.Color = [255 0 0];
            most_aligned_sam.Color = [255 0 0];
            pre = '<HTML><FONT color="';
            post = '</FONT></HTML>';
            listboxStr_ref = cell(numel( most_aligned_ref ),1);
            listboxStr_sam = cell(numel( most_aligned_sam ),1);
            
            str_ref = [pre reshape( dec2hex(most_aligned_ref.Color, 2)',1, 6) '">' most_aligned_ref.name post];
            listboxStr_ref = str_ref;
            Show_most_aligned_ref = {listboxStr_ref};
            str_sam = [pre reshape( dec2hex(most_aligned_sam.Color, 2)',1, 6) '">' most_aligned_sam.name post];
            listboxStr_sam = str_sam;
            Show_most_aligned_sam = {listboxStr_sam};
            
        else
            Show_most_aligned_ref = int2nt(mod_ref(n));
            Show_most_aligned_ref = cellstr(Show_most_aligned_ref);
            Show_most_aligned_sam = int2nt(most_pos_sam(n));
            Show_most_aligned_sam = cellstr(Show_most_aligned_sam);
        end
        
    Result_highest_quality(1,n) = Show_most_aligned_ref;
    Result_highest_quality(2,n) = Show_most_aligned_sam;
    
    end    
end

% Displaying qualities - Do not Edit!
fprintf('The qualities '); cprintf('_red', 'in the range'); fprintf(' are'); disp(quality_family(X));

fprintf('The cells '); cprintf('_red', '"Result_in_range"'); fprintf(' showing the possible alignment matching in the range of quality you had set. \n\n\n\n');

fprintf('The '); cprintf('_red', 'highest'); fprintf(' quality is'); disp(highest_q);

fprintf('The cells '); cprintf('_red', '"Result_highest_quality"'); fprintf(' showing the most possible alignment matching between refenrence gene and your sample gene. \n\n');