function rna_sequence_comparison
%searches for homologous subsequences in rna sequences

full_match=2;   %identical nt
half_match=1;   %conserved pyrimidine or purine
mis_match=0;    %purine to pyrimidine or vice-versa

nt_font_size=6; 

%specify sequence length and randomly generate sequences
seq_one_length=90;
seq_two_length=110;
sequence_one=rand(1,seq_one_length);
sequence_one=5*sequence_one;
sequence_one=ceil(sequence_one);
sequence_two=rand(seq_two_length,1);
sequence_two=4*sequence_two;
sequence_two=ceil(sequence_two);
%convert the random sequences into letters
sequence_one_alpha=num2str(sequence_one);
sequence_one_alpha=strrep(sequence_one_alpha,'1','a');
sequence_one_alpha=strrep(sequence_one_alpha,'2','c');
sequence_one_alpha=strrep(sequence_one_alpha,'3','g');
sequence_one_alpha=strrep(sequence_one_alpha,'4','u');
% sequence_one_alpha=strrep(sequence_one_alpha,'5','u');
% sequence_one_alpha=strrep(sequence_one_alpha,'6','u')
sequence_two_alpha=num2str(sequence_two');
sequence_two_alpha=strrep(sequence_two_alpha,'1','a');
sequence_two_alpha=strrep(sequence_two_alpha,'2','c');
sequence_two_alpha=strrep(sequence_two_alpha,'3','g');
sequence_two_alpha=strrep(sequence_two_alpha,'4','u');
% sequence_two_alpha=strrep(sequence_two_alpha,'5','u');
% sequence_two_alpha=strrep(sequence_two_alpha,'6','u')

%convert nucleotide sequences into appropriate numbers for pairwise sums
sequence_one=strrep(sequence_one_alpha,'a','1');
sequence_one=strrep(sequence_one,'u','2');
sequence_one=strrep(sequence_one,'g','4');
sequence_one=strrep(sequence_one,'c','8');
sequence_one=str2num(sequence_one);
sequence_two=strrep(sequence_two_alpha,'a','1');
sequence_two=strrep(sequence_two,'u','2');
sequence_two=strrep(sequence_two,'g','4');
sequence_two=strrep(sequence_two,'c','8');
sequence_two=str2num(sequence_two);

%make a matrix of the pairwise sums of each nt in the first sequence and each nt in the second sequence
seq_seq_matrix=[];
for i=1:length(sequence_one)
    for j=1:length(sequence_two)
        seq_seq_matrix(i,j)=sequence_one(i)+sequence_two(j);
    end
end
%size(seq_seq_matrix)

% use find and replace (using the table of which sums represent which pairs) in the seq_seq_matrix to change to the value for each pair's degree of matching
    %it is also possible to weight certain combinations more, eg cg correspondence is worth more than at correspondence
seq_seq_matrix(find(seq_seq_matrix==2))=full_match;
seq_seq_matrix(find(seq_seq_matrix==3))=half_match;
seq_seq_matrix(find(seq_seq_matrix==4))=full_match;
seq_seq_matrix(find(seq_seq_matrix==5))=mis_match;
seq_seq_matrix(find(seq_seq_matrix==6))=mis_match;
seq_seq_matrix(find(seq_seq_matrix==8))=full_match;
seq_seq_matrix(find(seq_seq_matrix==9))=mis_match;
seq_seq_matrix(find(seq_seq_matrix==10))=mis_match;
seq_seq_matrix(find(seq_seq_matrix==12))=half_match;
seq_seq_matrix(find(seq_seq_matrix==16))=full_match;

% starting in the first slot find each possible path (down and to the right are the valid directions) -- path length must be greater than one
% within the path find all groups of consecutive positive or consecutive negative numbers
% sum each group and (to reward longer matches or punish longer unmatches) raise it to the power of (one+ (the group's length)/ten) [again, these values need to be tuned]
% sum the total value for that path
% display the top n matches, showing the two matching subsequences, bolding the full matches and demibolding the half matches.  display the match score next to each match.

colormap(flipud(hsv(3)))  %3 values, 1 for each match condition -- flipping the color matrix sets red for the full match (when using negative value for mismatch will need to scale this)
imagesc(seq_seq_matrix)
hold on

text(repmat(1,1,length(find(sequence_one==1))),find(sequence_one==1),'a','fontsize',nt_font_size)
text(repmat(1,1,length(find(sequence_one==2))),find(sequence_one==2),'u','fontsize',nt_font_size)
text(repmat(1,1,length(find(sequence_one==4))),find(sequence_one==4),'g','fontsize',nt_font_size)
text(repmat(1,1,length(find(sequence_one==8))),find(sequence_one==8),'c','fontsize',nt_font_size)

text(find(sequence_two==1),repmat(length(sequence_one),1,length(find(sequence_two==1))),'a','fontsize',nt_font_size)
text(find(sequence_two==2),repmat(length(sequence_one),1,length(find(sequence_two==2))),'u','fontsize',nt_font_size)
text(find(sequence_two==4),repmat(length(sequence_one),1,length(find(sequence_two==4))),'g','fontsize',nt_font_size)
text(find(sequence_two==8),repmat(length(sequence_one),1,length(find(sequence_two==8))),'c','fontsize',nt_font_size)

% xlabel(sequence_two_alpha,'fontsize',8)  %ideal is to scale font size with string length, or draw the letter for 
% ylabel(sequence_one_alpha,'fontsize',6,'rotation',270) %rotate this label to have the letter corresponding to the match values it represents
set(gca,'xticklabel',[],'yticklabel',[])