%% Sequence Alignment, BLAST Search, and MSA
%% Download
%

%% Download Sequences - Myglobin 
mb_human = getgenbank ('BC014547', 'SequenceOnly', true);
mb_orangutan = getgenbank ('NM_001132084', 'SequenceOnly', true);
mb_mouse = getgenbank('BC025172', 'SequenceOnly', true);

%
%% Sequence Alignment - Global Sequence Alignment

% Aligning human and orangutan sequences for myoglobin
[s_ho, a_ho] = nwalign(mb_human, mb_orangutan, ...
                    'Alphabet', 'NT', 'ShowScore',true)

%

% Aligning human and mouse sequences for myoglobin
[s_hm, a_hm] = nwalign(mb_human, mb_mouse, ...
                    'Alphabet', 'NT', 'ShowScore',true)

%
%% Sequence Alignment - Local Sequence Alignment

% BRCA1 - breast cancer
seq_brc1a = getgenbank("U57574", "SequenceOnly", true);
% Immunoglobulin heavy variable
seq_ighv1 = getgenbank("NG_001019", "SequenceOnly", true);

sprintf('Sequence Lengths BRCA1: %d, and IGHV1: %d', ...
                 length(seq_brc1a), length(seq_ighv1))
[s_lcl, a_lcl] = swalign (seq_brc1a, seq_ighv1, 'Alphabet', 'NT');
disp(a_lcl);
sprintf('Alignment Score: %f5.2, length: %d', s_lcl, size(a_lcl, 2))
seqalignviewer(a_lcl)

%

% Aligning FOXO1 and FOXO3 Genes (Human)
% Comparing transcription factor binding sites 
% within promoter regions of related human genes

foxo_1 = getgenbank('NM_002015', 'SequenceOnly', true);
foxo_3 = getgenbank('NM_001455', 'SequenceOnly', true);

% Local alignment using swalign
[Score, Alignment] = swalign(foxo_1, foxo_3, 'Alphabet', 'NT');
disp(Score);
disp(Alignment);

%
%% BLAST Searches

S = getpdb('1CIV');
[RID1,ROTE] = blastncbi(S,'blastp','expect',1e-10);
report1 = getblast(RID1,'WaitTime',ROTE,'ToFile','1CIV_report.xml')

%

% Accession number of the hits and High Scoring Pairwise Segments
report1.Hits.Accession

%

[report1.Hits.Hsps]

%
%% Multiple Sequence Alignment

seq = {'ATTTACGCCT', 'TTAAGCCAT', 'TTAATTAACC', ...
                        'ATTTTCCGGA', 'AATTTACCGCCT'};
dist = seqpdist (seq, 'Alphabet', 'nt');

%

links = linkage(dist);
dendrogram (links);

%

seqs = {'AATTCCCGG' 'ATTCCGG' 'TTCCGG' 'AATTCTGG' 'AATTCC'};
ma = multialign(seqs) 
seqalignviewer(ma)

%
%% Consensus Sequence Determination

consensus = seqconsensus(ma, "Alphabet", "NT", "Gaps", "All")
% Remove gaps from the consensus
disp(consensus);
disp("Remove indels ...");
consensus = strrep(consensus, "-", "");
disp(consensus);
sprintf("Length of consensus is %d", strlength(consensus))