function gene_finder(name_of_file)
% GENE_FINDER(name_of_file) Run the program "Gene_finder"
%   name_of_file: name of FASTA file, e.g. 'name_of_file.fasta'

[header,seq] = readfasta(name_of_file);
sek = leading_and_complement(seq);

for j = 1:length(sek)
    seq = sek{j,1};
    
    for frame_offset = 1:3
        ends = find_ends_in_frame(seq, frame_offset);
        genes = find_genes_in_frame(seq, ends, frame_offset);
        
        genes_str = genes_to_string(genes, seq);
        
        freq_map = codon_frequency_map(seq, frame_offset);
        
        clear score;
        for i = 1:length(genes_str)
            gene_str = genes_str{i};
            
            scoring_windows = get_scoring_windows(gene_str);
            score(i) = gene_score(gene_str, scoring_windows, freq_map);
        end
        
        r{frame_offset} = apply_scoring(genes, score);
    end
    
    results{j} = r;
    
end
table_results = save_as_table(results, seq);

function [header, seq] = readfasta(filename)
% [HEADER, SEQ] = READFASTA(FILENAME)   Open and load sequence from FASTA
%                                       file
%   filename: name of file that contains the sequence
%   header: header of sequence
%   seq: sequence from file

data=importdata(filename);

header=data.Header;
seq=data.Sequence;

function sek = leading_and_complement(seq)
% SEK = LEADING_AND_COMPLEMENT(SEQ) Create a variable that contains
%                                   leading strand and its reverse complement
%   seq: leading strand
%   sek: variable of a type cell that contains both strands

sek{1,1} = seq;
sek{2,1} = seqrcomplement(seq);

function ends = find_ends_in_frame(seq, frame_offset)
% ENDS = FIND_END_IN_FRAME(SEQ, FRAME_OFFSET)   Returns vector of first
%                                               positions of stop codons
%   seq: sequence we search through
%   frame_offset: reading frame
%   ends: positions of first base in each stop codon

if (frame_offset <= 0 || frame_offset > 3)
    error('Incorrect reading frame offset');
end

seq_l = length(seq);
ends = 0;
j = 1;

for i = frame_offset:3:seq_l;
    if seq_l - i < 2
        break;
    end
    
    s = seq(i : i + 2);
    if strcmp(s, 'TAA') || strcmp(s, 'TAG') || strcmp(s, 'TGA')
        ends(j) = i;
        j = j + 1;
    end
end


function gene = find_gene_with_start(seq, start, pos_prev_end, pos_end)
% GENE = FIND_GENE_WITH_START(SEQ, START, POS_PREV_END, POS_END) Finds position of a specific start codon between two stop codons 
%   seq: sequence in which we look for the gene
%   start: start codon that we are looking for
%   pos_prev_end: position of previous stop codon
%   pos_end: position of stop codon

if strcmp(start, 'ATG') == 0 && strcmp(start, 'GTG') == 0 && strcmp(start, 'TTG') == 0
    error('Invalid start codon');
end

end_cur = seq(pos_end : pos_end + 2);
end_prev = seq(pos_prev_end : pos_prev_end + 2);

if strcmp(end_cur,'TAA') == 0 && strcmp(end_cur,'TAG') == 0 && strcmp(end_cur, 'TGA') == 0
    error('Invalid stop codon');
end

if pos_prev_end ~= 1 && strcmp(end_prev, 'TAA') == 0 && strcmp(end_prev, 'TAG') == 0 && strcmp(end_prev, 'TGA') == 0
    error('Invalid previous stop codon');
end

starts = strfind(seq(pos_prev_end : pos_end), start);

if isempty(starts)
    gene(1) = -1;
    gene(2) = -1;
    return;
end

gene_start = -1;
gene_end = -1;

for i = 1:length(starts)
    if mod(length((pos_prev_end + starts(i) - 1):(pos_end + 2)), 3) == 0
        gene_start = pos_prev_end + starts(i) - 1;
        gene_end = pos_end + 2;
        break
    end
end

gene(1) = gene_start;
gene(2) = gene_end;


function genes = find_genes_in_frame(seq, ends, frame_offset)
% GENES = FIND_GENES_IN_FRAME(SEQ, END, FRAME_OFFSET) Tries to find start of a gene for each given stop codon
%   seq: DNA strand
%   ends: positions of stop codons
%   frame_offset: reading frame
%   genes: list of gene positions along with information about reading frame

i = length(ends);
j = 1;
genes = [0 0 0];

while i >= 1
    if i ~= 1
        end_prev = ends(i - 1);
    else
        end_prev = 1;
    end
    
    gene = find_gene_with_start(seq, 'ATG', end_prev, ends(i));
    if (gene(1) == -1 && gene(2) == -1)
        gene = find_gene_with_start(seq, 'GTG', end_prev, ends(i));
        if (gene(1) == -1 && gene(2) == -1)
            gene = find_gene_with_start(seq, 'TTG', end_prev, ends(i));
        end
    end
    
    if (gene(1) ~= -1 && gene(2) ~= -1) && gene_length_ok(gene, 200)
        genes(j,:) = [gene, frame_offset];
        j = j + 1;
    end
    i = i - 1;
end
genes = flipud(genes);


function length_ok = gene_length_ok(gene, min_length)
% LENGTH_OK = GENE_LENGTH_OK(GENE, MIN_LENGTH)  Determine length of gene
%                                               and if gene is shorter than
%                                               minimum length it's deleted
%   gene: position of potencial genes
%   min_length: minimum length of gene
%   length_ok: return 1 if gene is long enough or return 0 if it's too short

l = gene(2) - gene(1) + 1;
length_ok = l > min_length;


function map = codon_frequency_map(seq, frame_offset)
% MAP = CODON_FREQUENCY_MAP(SEQ, FRAME_OFFSET) Creates mapping between codon name and its frequency score 
%   seq: sequence
%   frame_offset: reading frame
%   map: map from codon name to its frequency score

Codons = codoncount(seq,'Frame',frame_offset);
codons_names = fieldnames(Codons);

for i = 1:length(codons_names)
    a = codons_names{i};
    frequencies(i) = (Codons.(a))/(length(seq)/3)*1000;
    
end

map = containers.Map(codons_names, frequencies);


function scoring_windows = get_scoring_windows(gene_str)
% SCORING_WINDOWS = GET_SCORING_WINDOWS(GENE_STR)   Creates positions of scoring windows
%   gene_str: potential genes represented as string
%   scoring_windows: matrix of start positions of scoring windows and their lengths

window_size = 48;
windows = [0,0];

len = length(gene_str);
i = 1;
j = 1;

while i <= len
    if i + window_size > len
        n = len - i;
    else
        n = window_size;
    end
    
    windows(j, 1) = i;
    windows(j, 2) = n;
    
    j = j + 1;
    
    if n < window_size
        break
    end
    
    i = i + window_size / 2;
end
scoring_windows = windows;


function genes_string = genes_to_string(genes, seq)
% GENES_STRING = GENES_TO_STRING(GENES, SEQ)    Converts each gene to string                                           
%   genes: positions of potential genes
%   seq: sequence
%   genes_string: list of genes represented as strings

for i = 1:length(genes)
    genes_string{i} = seq(genes(i,1):genes(i,2));
end


function score = gene_score(gene_str, windows, freq_map)
% SCORE = GENE_SCORE(GENE_STR, WINDOWS, FREQ_MAP)   Determine score for a gene
%   gene_str: gene represented as a string
%   windows: matrix of start positions of scoring windows and their lengths
%   freq_map: map of each codon and its frequency
%   score: gene's score

s = ones(1, length(windows));

for i = 1:length(windows)
    for j = windows(i,1):3:(windows(i,1) + windows(i,2) - 1)
        codon = gene_str(j:j + 2);
        codon_score = freq_map(codon);
        
        s(i) = s(i) * codon_score;
    end
end

score = sum(s)/length(s);


function result = apply_scoring(genes, score)
% RESULT = APPLY_SCORING(GENES, SCORE) Decides whether potential gene is considered a real gene.
%   genes: positions of potential genes
%   score: gene's score
%   result: list of positions of final genes

j = 1;
m = mean(score);

if (length(genes) ~= length(score))
    error('Array must have same length')
end

for i = 1:length(score)
    if score(i) > m
        result(j,:) = genes(i,:);
        j = j + 1;
    end
end

function table_results = save_as_table(results, seq)
% TABLE_RESULTS = SAVE_AS_TABLE(RESULTS, SEQ)   Converts results to a table and saves it
%   results: list of positions of final genes
%   seq: sequence
%   table_results: list of gene's positions in table

final_leading = [results{1,1}{1,1};results{1,1}{1,2};results{1,1}{1,3}];
final_reversecomplement = [results{1,2}{1,1};results{1,2}{1,2};results{1,2}{1,3}];

final_rev = [0 0 0];

delka = length(seq);

for i = 1:length(final_reversecomplement)
    
    final_rev (i,1) = delka + 1 - final_reversecomplement(i,2);
    final_rev (i,2) = delka + 1 - final_reversecomplement(i,1);
    final_rev (i,3) = final_reversecomplement(i,3)*-1;
    
end

final = sortrows([final_leading;final_rev]);

for i = 1:length(final)
    final(i,4) = length(final(i,1):final(i,2));
end

gene_start = final(:,1);
gene_end = final (:,2);
reading_frame = final (:,3);
gene_length = final(:,4);
table_results = table(gene_start,gene_end,reading_frame,gene_length);
write(table_results,'results_of_gene_finder.xls');