%% Introduction to Matlab
% Table of Contents
% Resources
%     Videos
%     Simple plotting code
% Basics
%     Data Types
%         Vectors and Matrices
%     Plotting
%         Plotting and labeling
%     String comparisons and find function
%     Regular Expressions
%         Lookahead Operator: non-overlapping and overlapping search
% Functions
% Containers
%     Cell Arrays
%     Structures
%     Dictionaries
%     Tables
%     Reading Tables from a CSV or Excel file
%         Datamatrix
% File Operations
%     textscan()
% Review Questions
% Appendix
% Function Definitions
%%
%
%% Resources
%% Videos
% Introduction to Live Scripting
%

clear all
close
clc

%
%% Simple plotting code
%

x = 1:0.1:10
y = x .^ 2
plot (x,y)

%
%% Basics
%% Data Types
% Integers, Real, Strings, Booleans
% By default MATLAB creates numbers to be double precision
% The type of data in a variable is reported by the class command
%

% integers - can be int8, int16, int32, and int64
x = int32(9);
% real numbers -- single (32) or double (64) precision 
% explicitly defined as single precision
y = single(5.3);
% double precision is default number type
z = 1.5e-7;
% string variable
s = "This is a string";
% boolean variable true or false
b = true;

fprintf("class of x is %s and its value is %d", class(x), x);
fprintf("class of y is %s", class(y));
fprintf("class of b is %s", class(b));

%
%% Vectors and Matrices
% Vectors and Matrices are created using the square brackets notation. Elements are separated by commas or spaces, and a semicolon is designated a new row. The function length(v) returns the length of a vector. The function size(m) returns a vector with the dimensions of the matrix. The length of a matrix is the size of its longest dimension.
%

v = 1:10
sprintf("The length of vector v is %d", length(v))

m = [1 2 4 4; 8 1 6 5]
sprintf("The size of matrix m");
disp(size(m));

sprintf("The length of matrix m is %d", length(m))

%
%% Plotting
% Create a vector x of values from 0 to  in steps of 0.1
% Create another vector y by taking sin(x) -- i.e., trigonometric sine function for each value in x
% Plot y vs. x
%

% create x
x = 0:0.1:2*pi;
% create y
y = sin(x);
% plot y vs x
plot(x, y);
grid on
title("Sine Function")
xlabel("x")
ylabel("sin(x)")

%
%% Plotting and labeling
% Multiple plots are shown. One of the points is also labelled.


figure()
x = 0:0.1:2*pi;
y=sin(x);
plot(x,y)
grid on
% plotting the second graph
hold on
z = exp(-x);
plot(x,z,'r:*');
axis ([ 0 2*pi 0 1])
title('2-D Plots');
xlabel('Time');
ylabel('Sin(t)');

% This is text for creating the legend
text(pi/3,sin(pi/3),'\leftarrow (sin \pi/3)')
legend('Sine Wave','Decaying Exponential');

%
%% String comparisons and find function
% strfind function is used for finding the location of a pattern within a string.
% Ask the user to enter a string. Note that the user input should be enclosed in quotes for it to be treated as a string. Alternatively, one may use the option "s" to force the input to be treated as a string.
% Find the location of the pattern "the" in the string.
%

str = input("Enter a string to search: ", "s")
str = lower(str);
pat = "the";
locs = strfind(str, pat);
sprintf("Locations of %s ", pat);
disp(locs)

%
% The function strcmp is used to compare strings. The function strcmpi performs this comparison without considering the case of the strings.
% Create a string variable eco  and set to "Eco-RI"
% Create array enzymes of string variables = {"Bam-HI", "Eco-RI", "Hind-III", "Eco-RI", "Eco-RII"}
% Find the indices of eco within enzymes
% Find the indices of enzymes that contain the substring "RI"
%

eco = "Eco-Ri";
enzymes = ["Bam-HI", "Eco-RI", "Hind-III", "Eco-RI", "Eco-RII"];
% Find the indixes of eco within enzymes
match = strcmpi(eco, enzymes)
indx = find(match);
sprintf("Found %s at locations ", eco)
disp(indx)

% Find the indices of enzymes that contain the substring "RI"

indx = strfind(enzymes, "RI")   % create a cell array of matches
h = @(x) ~isempty(x);           % function handle
m = cellfun(h, indx)            % logical match array 
indx = find(m)                  % index of matches

% can combine all steps into a single command

fprintf ("executed as a single command \n");
find(cellfun(@(x) ~isempty(x), strfind(enzymes, "RI")))

%
%% Regular Expressions
% Regular expressions provide a powerful mechanism and syntax for specifying patterns that need to be searched. Regular expressions use special characters, such as "\s" for space characters, "\w" for word characters that include alphanumeric and underscore characters, and "\d" for any digits. A single character from a range of characters may be chosen by including the range in "[]".  For example, [aeiou] will match a vowel, and [^aeiou] will match any character except a vowel. A range of characters to choose from may also be specified, such as [A-Za-z] is a valid specification to choose an alphabetical character.
% Quantifiers are further used to enhance the number of times a token is matched. The quantifier * matches a token 0 to infinite times, + matches 1 to infinite times,? matches 0 or 1 times. Quantifier ranges may also supported with {m,n}, which matches the desired token m through n times.
% Some examples further illustrate the use of regular expressions:
% The regular expression, [aeiou]+ looks for the occurrence of one or more vowels
% The regular expression \s\w+y\s matches words surrounded by spaces and ending in a y
% The regular expression \w+ain[^\w] matches words that end in "ain"
% The regular expression [CG](TA){1,3} matches patterns where a C or G is followed by 1,2 or 3 occurrences of string TA.
%

str = "The rain in Spain gently falls in plain.";
pat = "[aeiou]+";
[s, m] = regexpi(str, pat, "start", "match")

%

pat = "\s\w+y\s";
[s, m] = regexpi(str, pat, "start", "match")

%

pat = "\w+ain[^\w]";
[s, m] = regexpi(str, pat, "start", "match")

%

str = "CCTATGTGTATACC";
pat = "[CG]TA{1,3}";
[s, m] = regexpi(str, pat, "start", "match")

%
%% Lookahead Operator: non-overlapping and overlapping search
% Consider the string "AAACCTAATAAA"
% Provide command for finding count of non-overlapping occurrence of "AA" or "AT" in the string.
% Provide command for finding the count of overlapping occurrence of the above pattern using the lookahead operator (?=pat), which matches pat but doesn’t consume it.
% As discussed above, the regular expression [AT] matches one character specified in the list of characters enclosed in square brackets.
%

seq = "AAACCTAATAAA";
freq = pattern_count("A[AT]", seq);
pat = "A[AT]";
loc = regexp(seq, pat);

fprintf ("Non-overlapping frequency count %-d", length(locs));

%

pat = "A(?=[AT])";
locs = regexp(seq, pat);
freq = pattern_count("A(?=[AT])", seq);
fprintf ("Overlapping frequency count %-d", length(locs));

%
%% Functions
% A function may be used to count patterns. Functions enable modularization and reusing the same logic.
% A function defined in the Appendix checks whether a set of 3 lengths forms a triangle. To check if a triangle is formed by three lengths, the triangle inequality must be satisfied. That is, the length of the largest side must be less than the sum of two other sides.
%

is_triangle(3, 4, 5)
is_triangle(3, 9, 5)

%
% Function pattern_count(defined at the end of live-script) computes the count of occurrences of a regular expression in a string.
% Alternatively, functions may be defined as separate ".m" files.
%

seq = "AAACCTAATAAA";
pat = "A[AT]";
freq = pattern_count(seq, pat);

fprintf ("Non-overlapping frequency count %-d", freq);

pat = "A(?=[AT])";
freq = pattern_count(seq, pat);
fprintf ("Overlapping frequency count %-d", freq);

%
%% Containers
%
%% Cell Arrays
% A cell array is a data type with indexed data containers called cells, where each cell can contain any type of data, including other cell arrays.
% Create an empty cell array with 3 rows and 2 columns
% The three rows of the cell array contain the following elements:
%             "eye", 3x3 identity matrix
%             "fibonacci,” array with elements: [1, 1, 2, 3, 5]
%             "magic,” a 3x3 magic square
% Display the cell array
%

mycell = cell(3,2);
mycell(1,1) = {"eye"};
mycell(1,2) = {eye(3)};

mycell(2,1) = {"fibonacci"};
mycell(2,2) = {[1, 1, 2, 3, 5]};

mycell(3,1) = {"magic"};
mycell(3,2) = {magic(3)};

disp ("Display elements of the cell")
celldisp(mycell)

disp ("Graphical Schematic of cell array contents")
cellplot(mycell)

%
%% Structures
% Create a structure array with patient records, with each patient's age, blood glucose, and heart rate;
%             John Smith, Age = 30, Blood Glucose = 109, Heart Rate = 80
%             Sally Snow, Age = 24, Blood Glucose = 99, Heart Rate = 68
%             Arun Seth, Age = 47, Blood Glucose=116, Heart Rate = 84
%             Sonia Bull, Age = 39, Blood Glucose=118, Heart Rate = 81
% Create a stacked bar graph with three series: Age, Glucose, and HR
%

% Create an array of structures

pats = struct();
pats(1).Name = "John Smith";
pats(1).Age = 30;
pats(1).Glucose = 109; 
pats(1).HR = 80;

pats(2).Name = "Sally Snow";
pats(2).Age = 24;
pats(2).Glucose = 99; 
pats(2).HR = 68;

pats(3).Name = "Arun Seth";
pats(3).Age = 47;
pats(3).Glucose = 116; 
pats(3).HR = 84;

pats(4).Name = "Sonia Bull";
pats(4).Age = 39;
pats(4).Glucose = 118; 
pats(4).HR = 81;

% extract ages, glucose, and heart rate

age = [pats.Age];
glucose = [pats.Glucose];
hr = [pats.HR];

figure();

res = [age' glucose' hr'];
names = categorical([pats.Name]);
bar(names, res);
grid;
legend({'Age', 'Glucose', 'HR'}, ...
    'Location', 'eastoutside');

%
%% Dictionaries
% Create a dictionary listing the locations of the vowels in the string: "The quick brown fox jumped over the lazy dog"
% Use the dictionary to list the match locations for the vowel "o"
%


vowels = ['a', 'e', 'i', 'o', 'u'];
vowels
vowels2 = ["a", "e", "i", "o", "u"];
vowels2

%

nums = [1 2 3 4 5]

%

for v = nums
    disp(v)
end

%

dict = dictionary();

str = "The quick brown fox jumped over the lazy dog";
vowels = ['a', 'e', 'i', 'o', 'u'];

for v = vowels
    indx = regexp(str, v, 'ignorecase');
    dict(v) = {indx};
end

disp("Locations of occurrence of 'o' ");
dict

%
%% Tables
% Create a cell array of names: "Adam,” "John,” "Sally"
% Create an array of scores: [81, 72, 90]
% Create a table of names and scores. The table's column names should be "Student", "Score"
% Add a new column to the table named "Grade.”
% Assign the grades to students based on logic:
%                 \text{Grade} =
% \left\{
%     \begin{array}{lr}
%         A, & \text{if score}  \geq 90\\
%         B, & \text{if }  80 \leq score  < 90 \\
%         C, & \text{if }  70 \leq score  < 80 \\
%         D, & \text{otherwise}
%     \end{array}
%

names = {"Adam", "John", "Sally"};
scores = [81, 72, 90];
class = table();
class.Student = names';
class.Score = scores';
disp(class);

% Loop through rows
for i = 1:size(class, 1)
    x = class.Score (i);
    if x >= 90
        class.Grade(i) = 'A';
    elseif x >= 80
        class.Grade(i) = 'B';
    elseif x >= 70
        class.Grade(i) = 'C';
    else
        class.Grade(i) = 'D';
    end
end

disp(class);

%
%% Reading Tables from a CSV or Excel file
% Read a file samples.csv into a table
% Compute the average strength value
% Add a new column named Positive, which is a Boolean (True or False) based on whether the Strength is greater than the average
% Save the revised table into a file samples_2.csv
%

samples = readtable ("samples.csv");
avg_strength = mean(samples.Strength);
samples.Postive = samples.Strength > avg_strength;
writetable(samples, "samples_2.csv");

%
%% Datamatrix
% A typical Datamatrix file is a CSV file with sample readings as its columns and genes along its rows. The first row is treated as VariableNames and may be used for accessing columns using a dot notation. For example, the first column is accessible as tbl.Gene.
% List the names of variables
% List the gene expression values for Sample_1 as a sub-table.
%

filename = "data_matrix.csv";
tbl = readtable (filename)
fprintf('matrix size %d %d', size(tbl));

fprintf("Names of the variables\n");
tbl.Properties.VariableNames

% Accessing sub-tables

disp('Accesing Gene Expression Values for Sample_1 ')
tbl(:, {'Gene', 'Sample_1'})

%
% Find the average expression value for Gene_3 across all the samples.
% The average needs to be computed using the values stored in numeric columns of the row for Gene_3. The numeric vector is a Boolean vector with true values for columns that are numeric, assuming that the first data row is representative of the table contents.
%

% check which of the variables are numeric
numeric = cellfun(@isnumeric, table2cell(tbl(1,:)))

%
% While table columns are accessible using names of the Variables, the relevant row for computing average expression value needs to be accessed by its row index. The row index is determined with strcmp().
%

% search the Gene column for "Gene_3"
row_index = strcmp(tbl.Gene, 'Gene_3')

%
% The average is next computed with numerical values in the row identified by row_index. Curly braces are used to access the contents of the table. Using parenthesis does not retrieve the table contents, but rather a sub-table wrapping the contents indexed.
%

avg = mean(tbl{row_index, numeric});
fprintf("Average for Gene 3 is %5.3f", avg);

%
% The first column is treated as RowNames if the 'ReadRowNames'option is set to true.
% Subtables and rows can then be accessed with "names" of the rows defined by the first column of the data matrix.
% Read values of Gene_4 for samples Sample_2 and Sample_5
% Find the value average expression of Gene_3 across all samples.
%

%% Reading and with ReadRowNames set
filename = "data_matrix.csv";
frm = readtable (filename, 'ReadRowNames', true);
disp('values of row Gene_4 in Sample_2 and Sample_5');
frm('Gene_4', {'Sample_2', 'Sample_5'})
% compute average expression for Gene_3

avg = mean(frm{"Gene_3", :});
sprintf("Average for Gene 3 is %5.3f", avg)

%
%% File Operations
% Text files may be read in line-by-line using the fgets function.
% For example, the following commands will read the file sample.csv and print out its contents. Note that the newline characters may often need to be removed from the end of the line. This is accomplished by using the strtrim() function. By enclosing the printed string within the delimiter "*" the code segment shows the line read from the file before and after trimming it.
%

fname = "sequence.fasta";
fid = fopen(fname);

% read past the header
txtline = fgets(fid);

seq = "";

%get the first line 
txtline = fgets(fid);
while ischar(txtline) 
    %fprintf ("*%s*", txtline);
    %fprintf ("%s\n", strtrim(txtline));
    seq = strcat(seq, strtrim(txtline));
    txtline = fgets(fid);
end
fclose(fid);

%

seq = read_fasta_seq("sequence.fasta");
disp(seq)

%
%% textscan()
% Using the textscan() function enables reading text data in a custom format from a character stream or a file stream. While the full capabilities of this function are listed here, textscan() allows the specification of custom delimiters, specialized treatment of empty and NaN values, as well as skipping comments and literal text.
% Generally, read data requires that the format of the expected data types be provided. textscan then reads the data into a cell array with the number of elements in a row being equal to the number of format specifiers and data being read getting stacked as columns within each cell.
% The following example used the textscan function to read numbers from a string. The first example shown reads only the first two floating point numbers from the string since the format specifier is used for the specified frequency of 2 times. When the frequency is not specified, the format specifier is used repeatedly, and all four numbers are read in. The format specifier determines the number of columns in the cell array returned. The elements are read into a cell array with three columns as in the last example.
%

str = "0.75 8.9 5.334 1e-3";
val = textscan(str, "%f", 2);
val{1}

val = textscan(str, "%f");
val{1}

val = textscan(str, "%f %f %f")
val{1}

%
% The textscan()  function can similarly read values from a file with the file pointer specified as the string or stream to be read in. The number of expected to be read in can be used to create a format string using the repmat construct which can create a string of format specifiers. With the proper number of elements in the format specifier, the elements are read as a data matrix.
%

%% Open File 
fid=fopen('data_matrix.csv');

numCols = 6;
numHeaders = 1;
% Define Formats
format_header = repmat('%s ', 1, numCols);
format_data = ['%s ', repmat('%f ', 1, numCols-1)];

% Read Data using defined formats
header = textscan(fid, format_header, 1);
data = textscan(fid, format_data);
% Close file
fclose(fid);

% sample names
disp("sample names"); 
[header{2:numCols}]

%gene names
disp ("gene names"); 
[data{:,1}]

%expression data
disp ("expression data");
[data{:, 2:numCols}]

disp ('Elements of first and last columns');
[first,~, ~, ~, ~, ~, last] = deal(data)

%
%% Review Questions
%
% Question 1: Plot the  function in the range of to by stepping through the  values in increments of 0.05.
%

hold off 
x = 0:0.05:2*pi
y = cos(x)
plot(x,y)

%
% Question 2: A string is an array of characters. Search and list the index locations of character "A" in the string TTATTATTATTA.
%

str = "TTATTATTATTA";
pat = "A";
s_find = strfind(str,pat)

%
% Question 3:
% Consider the following array of strings:
% ["AATAA", "CCTG", "AAA", "AGGCT", "AAAT"]
% (a) Find the indices of elements that are equal to the string "AA"
% (b) Find the indices of elements that include the pattern "AA"
%

str = ["AATAA", "CCTG", "AAA", "AGGCT", "AAAT"];
pat = "AA";
%(a)
find(strcmp(str, pat))

%

%(b)
str = ["AATAA", "CCTG", "AAA", "AGGCT", "AAAT"];
pat = "AA";
indx = strfind(str, pat);
find(cellfun(@(x) ~isempty(x), indx))

%
% Question 4:
% Consider the string:
% "She sells sea shells on the sea shore."
% Use regular expressions to find:
% (a) all occurrences where a vowel is followed by another vowel
% (b) A vowel is followed by exactly two occurrences of the letter "l" (el).
% (c) Words that end in the letter "e".
%


str = "She sells sea shells on the sea shore.";
%(a)
pat = "[aeiou]{2,2}"
[s, m] = regexpi(str, pat, 'start', 'match')

%(b) 
pat = "[aieou](l){2,2}"
[s, m] = regexpi(str, pat, 'start', 'match')

%(c)
pat = "\w+e[^\w]";
[s,m] = regexpi(str, pat, "start", "match")

%
% Question 5: Find the overlapping and non-overlapping frequency counts for the pattern "CC" or "CG" in the sequence CCCGCCGGCGCCCC
%

seq = "CCCGCCGGCGCCCC"
% Non overlapping patterns CC or CG
re = "C[CG]"
regexpi(seq, re)

% overlapping hits
re= "C(?=[CG])"
regexpi(seq, re)

%
% Question 6: A stop codon in a gene sequence may be either "TAG," "TGA," or a "TAA." Visit the GenBank and download the DNA sequence having the accession number of M10051, saving it as a FASTA formatted record. List the locations of a stop codon in the sequence saved. (Note that the first line in the saved file should be ignored).
%

seq = read_fasta_seq("sequence.fasta");

re = "T(AG|GA|AA)";
[s, m] = regexpi(seq, re, "start", "match")

%
% Question 7: A sequence is AT-rich when the total count of nucleotides (A+T) is more than 50%. Write a function that returns the percentage of nucleotides found in the input sequence. The prototype of the function:
% [frac] = function at_fraction (seq)
% Test your function by calling it with the following input sequence having a (A+T) fraction of 60%.
% "AATTGGCCAATTCCCCAATA"
%
% Question 8: In the patient array of structures, add two new patients in Structures
%             Ian Brown, Age = 51, Blood Glucose = 105, Heart Rate = 78
%             Sam Singh, Age = 36, Blood Glucose = 95, Heart Rate = 72
% Plot the bar graphs. Also, compute the mean value for blood glucose for the entire group.
%



%
% Question 9: With reference to the table created by reading the file data_matrix.csv in section Datamatrix
% create a new column named Ratio_1_2 with values equal to the ratio of gene expressions equal to the value for Sample_1 divided by the value for Sample_2.
% create another column named AvgExpr with a value equal to the average value of gene expression across all samples
%



%
% Question 10: For the cell array example in the section Cell Arrays
% replace the 3x3 magic square with a 4x4 magic square.
% add a third column with a short description for the rows: (1) "A 3x3 identity matrix", (2) "first five elements of the Fibonacci series", (3) a 4x4 magic square
% add fourth and fifth rows with contents:
%                 "zeros", zeros(3), "a 3x3 zero matrix"
%                 "ones", ones(4, 5), "a 4 by 5 matrix of ones"
%

mycell = cell(1,1);
mycell(1,1) = {"eye"};
mycell(1,2) = {eye(3)};

mycell(2,1) = {"fibonacci"};
mycell(2,2) = {[1, 1, 2, 3, 5]};

mycell(3,1) = {"magic"};
mycell(3,2) = {magic(3)};

%replaing element {3,2}
mycell{3,2} = {magic(4)};
mycell(1,3) = {"A 3x3 identity matrix"};
mycell{2,3} = {"first five elements of the Fibonacci series"};
mycell{3,3} = {"a 4x4 magic square"};

mycell{4,1} = {"zeros"};
mycell{4,2} = {zeros(3)};
mycell{4,3} = {"a 3x3 zero matrix"};

mycell{5,1} = {"ones"};
mycell{5,2} = {ones(4,5)};
mycell{5,3} = {"a 4 by 5 matrix of ones"};

disp ("Display elements of the cell")
celldisp(mycell)

disp ("Graphical Schematic of cell array contents")
cellplot(mycell)

%
% Thiss is a test for somee thing with spelling errors..
%
%% Appendix
% Functions can either be defined as .m files or appear at the end of a Live Script.
%% Function Definitions
%
% Reads a fasta file
%

function [seq] = read_fasta_seq(fname)
    fid = fopen(fname);

    % read past the header
    fgets(fid);

    seq = "";

    %get the first line 
    txtline = fgets(fid);
    while ischar(txtline) 
        seq = strcat(seq, strtrim(txtline));
        txtline = fgets(fid);
    end
    fclose(fid);
end

% Computes the frequency of a pattern in a sequence. Pattern and sequence is passed in as parameters. pattern is specified as a regular expression.


function [freq_count, locs] = pattern_count (sequence, pattern)
    locs = regexp(sequence, pattern);
    freq_count = length(locs);
end

% Determines if the three lengths a, b, c passed in as parameters could form a triangle.

function str_tri_test = is_triangle(a,b,c)
% istriangle determine if a,b,c form a triangle
    v = [a b c];
    vsort = sort(v);
    
    if vsort(1) + vsort(2) > vsort(3)
    s1 = ['The sides ', num2str(v),' form a triangle.'];
    % Checking the type of triangle formed
    if (a==b) & (b==c) % all sides equal
        s2 =' The triangle is an equilateral triangle.';
    elseif (a==b) | (b==c) | (a==c)
        s2 = ' The triangle is an isosceles triangle.';
    else
        s2 = ' The triangle is a scalene triangle.';
    end    
    str_tri_test = strcat(s1,s2);
    else
    str_tri_test = ['The sides ', num2str(v), ...
				       ' do not form a triangle.'];
    end
end