%EPS 109 Computer Simulations in the Earth Sciences
%December 2 2015
%Juneyoung Jeong

%Chaos Game Representation (CGR) of a Genome

%%Introduction: This is a method of visualizing an entire genome structure,
%%while maintaining information on the order of nucleotides. 
%%This idea was first introduced in a paper I found (Jeffrey, H., 1990) when choosing a topic. 

%%Basic Rules to Play:
%%1. Start at the center of a polygon
%%2. Randomly choose one corner (i.e. roll a dice)
%%3. Plot midpoint between the corner and where you currently are
%%4. Repeat steps 2-3 several thousand times

%%To summarize the results, the Chaos Game creates interesting shapes when 
%%the polygon has 3, 5, 6, or 7 vertices. These images/attractors 
%%exhibit periodicity, but other polygons are simply filled in, and are
%%chaotic.

%%This chaotic behavior also occurs with quadrilaterals. CGR applies this
%%knowledge to model gene sequences to suggest some sort of structure/order
%%in the sequencing. Instead of rolling a die, each corner represents a 
%%nucleotide, and the Chaos Game is played as if the DNA sequence is a list
%%of the resulting dice rolls.

%%This is illustrated with sparse gaps in the box.
%%Jeffrey observed that some of these gaps are characteristic in shape and size for similar
%%genomes (i.e. vertebrates have a distinct pattern, supposedly).



%% Initial Conditions
clear all
clf 

global seq corner_g corner_t corner_a corner_c

%unit square, upper-right corner at (1,1)
%upper-right corner, and clockwise, GTAC
corner_g=[1,1];
corner_t=[1,0];
corner_a=[0,0];
corner_c=[0,1];

%initial point in the center of the square
center=[0.5,0.5];

f=0; %image counter

%% Given Sequence; simply replace filename to read new sequence
seq=fileread('HUMHBB Sequence.txt');


%% Playing the Chaos Game, with our DNA dice
P = center; %initial point is center of unit square
for i=1:length(seq)
%     P = newPos(i,P);
    P = chaosPos(i,P); %for random choice, i.e. w/o DNA
    xx(i)=P(1);
    yy(i)=P(2);
    if mod(i,100)==0 %based on genome length. change accordingly
    plot(xx,yy,'.', 'MarkerSize', 4);
    pause(0.001);
%     f=f+1;
%     saveas(gcf, ['CGRform_',sprintf('%03d',f),'.png'])
    end
end
plot(xx,yy,'.', 'MarkerSize', 4);
% saveas(gcf, ['CGRform_',sprintf('%03d',f),'.png'])


%To observe fractal formation, gradually reduce MarkerSize
f=0;
figure
for j=1:20
    dotsize=8-j/5;
    plot(xx,yy,'.', 'MarkerSize', dotsize);
    axis([0 1 0 1]);
%     f=f+1;
%     saveas(gcf, ['CGRmarker_',sprintf('%03d',f),'.png'])
    pause(1);
end
