function [empty_genes, missing_gene_in_sequence, variable_parts, variable_parts_structure] = find_variable_parts(number_of_sequences, number_of_melt)
% FIND_VARIABLE_PARTS(number_of_sequences, number_of_melt) Run the program "Find_variable_parts"
%   number_of_sequences: number of genomes in analyzed files
%   number_of_melt: number of melt-types for analyzed genomes

empty_genes = {};
missing_gene_in_sequence = {};

[genes, number_of_genes] = fasta_from_current_dir;
variable_parts = {};
variable_parts_structure = struct;

melt{:,1} = {'S2';'S5';'S7';'S9';'S11';'S12'};
melt{:,2} = {'S1';'S3';'S4';'S6';'S8';'S10'};
melt{:,3} = {'S13';'S14';'S15';'S16';'S17';'S18'};
melt{:,4} = {'S19';'S20';'S21';'S22';'S23';'S24'};

for gene_number = 1:number_of_genes
    
    [sequence, header, number_of_seq,empty_genes] = load_sequences(genes,gene_number,empty_genes);
    [gene_in_all_sequences,missing_gene_in_sequence] = is_gene_in_all_sequences(number_of_sequences,number_of_seq,genes,gene_number,missing_gene_in_sequence);
    missing_gene_in_sequence = setdiff(missing_gene_in_sequence,empty_genes);
    
    if not(gene_in_all_sequences)
        continue
    end
    
    alignment = multiple_alignment_of_sequence(genes, gene_number);
    
    if not(is_cut_effective(alignment))
        continue
    end
    
    if neccessary_to_cut_sequences(sequence, number_of_seq)
        alignments = sequence_cut(alignment);
    else
        name = ['alignment_' num2str(1)];
        alignments.(name) = alignment;
    end
    
    alignments = alignment_length_for_alignment(alignments);
    
    if isempty(alignments)
        continue
    end
    
    for l = 1:numel(fieldnames(alignments))
        
        if not(count_procentual_changes_for_alignment(alignments, l))
            continue
        end
        
        name = ['alignment_' num2str(l)];
        alignment_index = count_alignment_index(alignments.(name));
        
        if isempty(alignment_index)
            continue
        end
        
        sequence_for_variability = sequence_for_analysis(alignment_index, alignments.(name), number_of_seq,header);
        distinguish_between_clusters = phylogenetic_analysis(sequence_for_variability, number_of_melt, header, number_of_seq, melt);
        
        for j = 1:length(distinguish_between_clusters)
            if distinguish_between_clusters(j)
                if melting_temperatures(sequence_for_variability, number_of_melt, header, number_of_seq, j, melt)
                    name = ['sequence_for_index_' num2str(j)];
                    variable_parts{end+1} = [genes(gene_number).name ' ' name];
                    index = strfind(genes(gene_number).name,'.')-1;
                    name_for_struct = [genes(gene_number).name(1:index) '_' name];
                    variable_parts_structure.(name_for_struct) = sequence_for_variability.(name);
                end
            end
        end
        
    end
end


function [genes, number_of_genes] = fasta_from_current_dir
% [GENES, NUMBER_OF_GENES] = FASTA_FROM_CURRENT_DIR  Load all FASTA files
%                                                    to structure
%   genes: structure contains information about FASTA files
%   number_of_genes: number of FASTA files

genes = dir('*.fasta');
number_of_genes = length(genes);

function [sequences, headers, number_of_seq,empty_genes] = load_sequences(genes,gene_number,empty_genes)
% [SEQUENCES, HEADERS, NUMBER_OF_SEQ, EMPTY_GENES] = LOAD_SEQUENCES(GENES,GENE_NUMBER,EMPTY_GENES)  Get sequences from FASTA files
%   genes: structure contains information about FASTA files
%   gene_number: number of FASTA files
%   empty_genes: predefined variable
%   sequences: sequences from file
%   headers: header of sequence
%   number_of_seq:number of sequences in file
%   empty_genes: genes from reference sequence which are not found in
%                genome of particular strains

if (genes(gene_number).bytes) ~= 0                                          
    [headers,sequences] = fastaread(genes(gene_number).name);
    number_of_seq = length(sequences);
else
    sequences = '';
    headers = '';
    number_of_seq = [];
    empty_genes{end+1} = genes(gene_number).name;
end

function [gene_in_all_sequences,missing_gene_in_sequence] = is_gene_in_all_sequences(number_of_sequences,number_of_seq,genes,gene_number,missing_gene_in_sequence)
% [GENE_IN_ALL_SEQUENCES, MISSING_GENE_IN_SEQUENCE] = IS_GENE_IN_ALL_SEQUENCES(NUMBER_OF_SEQUENCES,NUMBER_OF_SEQ,GENES,GENE_NUMBER,MISSING_GENE_IN_SEQUENCE) Check if gene is in all analyzed genomes
%   number_of_sequences: number of genomes in analyzed files
%   number_of_seq: number of sequences in file
%   genes: structure contains information about FASTA files
%   gene_number: number of FASTA files
%   missing_gene_in_sequence: predefined variable
%   gene_in_all_sequences: return 1 if gene is in all sequences or return 0
%   missing_gene_in_sequence: list of genes which are not in all sequences

if number_of_sequences == number_of_seq
    gene_in_all_sequences = 1;
else
    gene_in_all_sequences = 0;
    missing_gene_in_sequence{end+1} = genes(gene_number).name;
end

function alignment = multiple_alignment_of_sequence(genes, gene_number)
% [ALIGNMENT] = MULTIPLE_ALIGNMENT_OF_SEQUENCES(GENES,GENE_NUMBER) Load gene alignment
%   genes: structure contains information about FASTA files
%   gene_number: number of FASTA files
%   alignment: multiple alignment of gene

index = strfind(genes(gene_number).name,'.')-1;
name = [genes(gene_number).name(1:index) '_alignment.aln'];
[~, alignment] = multialignread(name);
alignment = char(alignment);

% alignment = multialign(sequence,'ScoringMatrix','NUC44');

function cut_alignment = neccessary_to_cut_sequences(sequence, number_of_seq)
% CUT_ALIGNMENT = NECCESSARY_TO_CUT_SEQUENCES(SEQUENCE, NUMBER_OF_SEQ) Check if is neccessary to cut sequences
%   sequence: sequences from file
%   number_of_seq: number of sequences in file
%   cut_alignment: return 1 if is neccessary to cut sequences or return 0

sequences_length = zeros(1,number_of_seq);

for i = 1:number_of_seq
    sequences_length(i) = length(sequence{i});
end

if length(unique(sequences_length)) == 1
    cut_alignment = 0;
else
    cut_alignment = 1;
end

function effective = is_cut_effective(alignment)
% EFFECTIVE = IS_CUT_EFFECTIVE(ALIGNMENT) Define effective cut of sequence as a part of genome that is well
%                                         defined (doesn't contain -) in all strains
%   alignment: multiple alignment of gene
%   effective: return 1 if is effective to cut sequences or return 0

defined = zeros(1,size(alignment,2));
effective = 0;

for i = 1:size(alignment,2)
    if sum((alignment(:,i)) == '-') == 0
        defined(i) = 1;        
    end
end

if sum(defined == 1) > 0
    effective = 1;
end

function alignment_cut = sequence_cut(alignment)
% ALIGNMENT_CUT = SEQUENCE_CUT(ALIGNMENT) Cut alignment of gene
%   alignment: multiple alignment of gene
%   alignment_cut: alignment of gene which not contains any '-'

defined = zeros(1,size(alignment,2));

for i = 1:size(alignment,2)
    if sum((alignment(:,i)) == '-') == 0
        defined(i) = 1;
    end
end

starts_of_defined_positions = find(diff([0 defined])==1);     
ends_of_defined_positions = [];

for i = starts_of_defined_positions    
    
    j = i;
    while j <= length(defined) && (defined(j) == 1) 
        j = j+1;
    end
    ends_of_defined_positions(end+1) = j-1; 
end

for j = 1:length(starts_of_defined_positions)
    
    name{j} = ['alignment_' num2str(j)];
    alignment_cut.(name{j}) = alignment(:,starts_of_defined_positions(j):ends_of_defined_positions(j));

end

function alignments = alignment_length_for_alignment(alignment)
% ALIGNMENTS = ALIGNMENT_LENGTH_FOR_ALIGNMENT(ALIGNMENT)  Determine length
%                                                         of alignment
%   alignment: multiple alignment of gene
%   alignments: length of alignment of gene

number_of_alignment = numel(fieldnames(alignment));
index = 1;

for i = 1:number_of_alignment
    name = ['alignment_' num2str(i)];
    align = alignment.(name);
    
    if alignment_length(align)
        names = ['alignment_' num2str(index)];
        alignments.(names) = align;
        index = index + 1;
    end
end

if exist('alignments') == 0
    alignments = [];
end

function long_enough = alignment_length(align)
% BOOL = ALIGNMENT_LENGTH(ALIGN) Check if alignment is longer than 70 bp.
%   align: alignment of sequences
%   long_enough: return 1 if is longer than 70 bp or return 0

if size(align,2) < 70
    long_enough = 0;
else
    long_enough = 1;
end

function variable_enough = count_procentual_changes_for_alignment(alignment, l)
% VARIABLE_ENOUGH = COUNT_PROCENTUAL_CHANGES_FOR_ALIGNMENT(ALIGNMENT, L) Count procentual of variable positions in alignment
%   alignment: multiple alignment of gene
%   l: index of alignment in file
%   variable_enough: return 1 if is contains more than 3% of variable positions or return 0

name = ['alignment_' num2str(l)];
align = alignment.(name);

change = 0;

for k = 1:size(align,2)
    
    if length(unique(align(:,k))) > 1
        change = change + 1;
    end
end

procentual_change = (change/size(align,2)) * 100;

if procentual_change >= 3
    variable_enough = 1;
else
    variable_enough = 0;
end

function alignment_index = count_alignment_index(alignment)
% ALIGNMENT_INDEX = COUNT_ALIGNMENT_INDEX(ALIGNMENT) Get index of variable
%                                                    parts in alignment
%   alignment: multiple alignment of gene
%   alignment_index: index of variable part in gene alignment

window = 100;
change = 0;
alignment_index = [];

if size(alignment,2) >= 70  && size(alignment,2)< 100 
    
    for i = 1:size(alignment,2)
        if length(unique(alignment(:,i))) > 1
            change = change + 1;
        end
    end
    
    if change >= 10
        alignment_index(end+1,:) = [1 size(alignment,2)];
    end
    
else
    
    for i = 1:size(alignment,2)-100+1
        aln = alignment(:,i:i+window-1);
        
        for k = 1:size(aln,2)
            if length(unique(aln(:,k))) > 1
                change = change + 1;
            end
        end
        
        if change >= 10
            position(i) = 1;
            change = 0;
        else
            position(i) = 0;
            change = 0;
        end
        
    end
    
    index = 1;
    
    if position(1) == 1
        window2(index,1) = 1;
    end
    
    for i = 2:length(position) - 1
        if position(i) == 1 && position(i-1) == 0
            window2(index,1) = i;
        end
        
        if position(i) == 1 && position(i+1) == 0
            window2(index,2) = i;
            index = index + 1;
        end
    end
    
    if exist('window2')
        
        if length(window2) == 1
            if (sum(position(window2:length(position)) == 1)) == (length(window2:length(position)))
                window2(1,2) = length(position);
            end
        else
            [number_of_rows,number_of_columns] = size(window2);
            if window2(number_of_rows,number_of_columns) == 0
                window2(number_of_rows,number_of_columns) = length(position);
            end
        end
        
        for i = 1:size(window2,1)
            alignment_index(end+1,:) = [window2(i,1) window2(i,2)+window-1];
        end
        
    end
    
    i = 1;
       
    while i < size(alignment_index,1)
        if alignment_index(i,2) > alignment_index(i+1,1)
            alignment_index(i,:) = [alignment_index(i,1) alignment_index(i+1,2)];
            alignment_index(i+1,:) = [];
        else
            i = i+1;
        end
        
    end
    
end

function sequence = sequence_for_analysis(alignment_index, alignment, number_of_seq, header)
% SEQUENCE = SEQUENCE_FOR_ANALYSIS(ALIGNMENT_INDEX, ALIGNMENT,
% NUMBER_OF_SEQ, HEADER) Get sequence from alignment.
%   alignment_index: index of variable part in gene alignment
%   alignment: multiple alignment of gene
%   number_of_seq: number of sequences in file
%   header: header of sequence
%   sequence: variable parts of sequences from alignment

sequence = struct;

for k = 1:size(alignment_index,1)
    
    for j = 1:number_of_seq
        
        name = ['sequence_for_index_' num2str(k)];
        name2_cell = header(j);
        name2 = name2_cell{1};
        
        sequence.(name).(name2) = alignment(j,alignment_index(k,1):alignment_index(k,2));
        
    end
end  

function distinguish_between_clusters = phylogenetic_analysis(sequence, number_of_melt, header, number_of_seq, melt)
% DISTINGUISH_BETWEEN_CLUSTERS = PHYLOGENETIC_ANALYSIS(SEQUENCE, NUMBER_OF_MELT, HEADER,
% NUMBER_OF_SEQ, MELT) Perform phylogenetic analysis for variable parts and
% check if it is possible to distinguish between individual melt-types.
%   sequence: variable parts of sequences from alignment
%   number_of_melt: number of melt-types for analyzed genomes
%   header: header of sequence
%   number_of_seq: number of sequences in file
%   melt: structure with header of sequences with belong to same melt-types
%   distinguish_between_clusters: return 1 if is possible to distinguish between individual melt-types or return 0

number_of_parts = numel(fieldnames(sequence));

for i = 1:number_of_parts
    
    name = ['sequence_for_index_' num2str(i)];
    seq = get_sequence_to_cell(sequence, name, number_of_seq, header);
    
    if N_in_sequences(seq)
        distinguish_between_clusters(i) = 0;
         continue
    end
    
    p_distance = seqpdist(seq,'Alphabet','NT');
    PhyloTree = seqlinkage(p_distance/2, 'average', header);
    clusters = cluster(PhyloTree,[],'criterion','maximum','maxclust',number_of_melt);
    
    unique_char = unique(clusters);
    number_of_occurence_in_cluster = zeros(1,length(unique_char));
    
    if  length(unique_char) == number_of_melt
        
        for k = 1:length(unique_char)
            number_of_occurence_in_cluster(k) = count(num2str(clusters'),num2str(unique_char(k)));
        end
        
        if length(unique(number_of_occurence_in_cluster)) == 1
            distinguish_between_clusters(i) = compare_cluster(header,clusters,melt);
        else
            distinguish_between_clusters(i) = 0;
        end
        
    else
        distinguish_between_clusters(i) = 0;
    end
end

function seq = get_sequence_to_cell(sequence, name, number_of_seq, header)
% SEQ = GET_SEQUENCE_TO_CELL(SEQUENCE, NAME, NUMBER_OF_SEQ, HEADER) Convert sequences to cell
%   sequence: variable parts of sequences from alignment
%   name: name of sequences from variable alignment
%   number_of_seq: number of sequences in file
%   header: header of sequence
%   seq: sequences in cell

for i = 1:number_of_seq
   
    name2 = header{i};
    seq{i} = sequence.(name).(name2);
    
end

function contains_N = N_in_sequences(seq)
% CONTAINS_N = N_IN_SEQUENCES(SEQ) Check if sequences contain base 'N'
%   seq: multiple alignment of gene
%   contains_N: return 1 if sequences contain 'N' or return 0

seq = upper(seq);
N_content = [];

for i = 1:length(seq)
   
   if sum(seq{i}=='N') >= 1
       N_content(end+1) = 1;
   else
       N_content(end+1) = 0;
   end
    
end

if sum(N_content) >= 1
    contains_N = 1;
else
    contains_N = 0;
end

function distinguish_between_clusters = melting_temperatures(sequence, number_of_melt, header, number_of_seq, j, melt)
% DISTINGUISH_BETWEEN_CLUSTERS = MELTING_TEMPERATURES(SEQUENCE, NUMBER_OF_MELT, HEADER, NUMBER_OF_SEQ, J, MELT) Perform calculation of melting temperatures for variable parts and
%                                                                                       check if it is possible to distinguish between individual melt-types.
%   sequence: variable parts of sequences from alignment
%   number_of_melt: number of melt-types for analyzed genomes
%   header: header of sequence
%   number_of_seq: number of sequences in file
%   j: index of variable part of alignment
%   melt: structure with header of sequences with belong to same melt-types
%   distinguish_between_clusters: return 1 if is possible to distinguish between individual melt-types or return 0

name = ['sequence_for_index_' num2str(j)];
seq = get_sequence_to_cell(sequence, name, number_of_seq, header);

seq = upper(seq);

load('nearest_neighbor_values.mat');

A = -0.0108;
R = 0.00199;
C = 0.0000005;
Na = 0.05;
H = 0;
S = 0;

temperature = zeros(1,24);

for j = 1:number_of_seq
    
    for i = 1:length(seq{j})-1
        
        row_index = ismember(T.nearest_neighbors,(seq{j}(i:i+1)));
        
        if sum(row_index == 1)
            
            H = H + T.delta_H(row_index);
            S = S + T.delta_S(row_index);
        else
            disp('unknown neighbour')
        end
        
    end
    temperature(1,j) = H/(A+S+R*log(C/4))-273.15+16.6*log10(Na);
    H = 0;
    S = 0;
end

tree = linkage(temperature','average');
clusters = cluster(tree,'maxclust',number_of_melt);
unique_char = unique(clusters);
number_of_occurence_in_cluster = zeros(1,length(unique_char));

if  length(unique_char) == number_of_melt
    
    for k = 1:length(unique_char)
        number_of_occurence_in_cluster(k) = count(num2str(clusters'),num2str(unique_char(k)));
    end
    
    if length(unique(number_of_occurence_in_cluster)) == 1
        distinguish_between_clusters = compare_cluster(header,clusters,melt);
    else
        distinguish_between_clusters = 0;
    end
    
else
    distinguish_between_clusters = 0;
    
end

function bool = compare_cluster(header,clusters,melt)
% BOOL = COMPARE_CLUSTER(HEADER,CLUSTERS,MELT) Compare cluster obtained from phylogenetic or melting temperature analysis with set one.
%   header: header of sequence
%   clusters: clusters contains single sequences 
%   melt: structure with header of sequences with belong to same melt-types
%   bool: return 1 if clusters are identical or return 0

header = header';

melt_from_cluster{:,1} = header(clusters == 1);
melt_from_cluster{:,2} = header(clusters == 2);
melt_from_cluster{:,3} = header(clusters == 3);
melt_from_cluster{:,4} = header(clusters == 4);

number_of_correct_clusters = [];

for i = 1:4
    for j = 1:4
        
        if isempty(setdiff(melt_from_cluster{1,i},melt{1,j})) && isempty(setdiff(melt{1,j},melt_from_cluster{1,i})) 
            number_of_correct_clusters(end+1) = 1;
        end

    end
end

if length(number_of_correct_clusters) == 4
    bool = 1;
else
    bool = 0;
end


