function do_Mus_VDJ_alignments(data_filename,start, skip)
%%%%%%%%%%%%%%
% This function reads the data file with sequence reads;
% produces, for each raw NT sequence, a 'minimal' list of plausible
% VDJ recombination events, described by gene choices,
% deletions, insertions, etc.
%
% By 'minimal', we mean that only the information necessary to
% generate a full list of plausible events is stored in the output
% files -- for instance, only the minimum number of deletions
% and insertions are stored. The expectation maximization algorithm
% reads these 'alignment' files and loops over all events (greater
% numbers of deletions and insertions) using this information.
%

%% Load Data
if ~exist('start','var')
    start = 1; % use every nth sequence.
end
if ~exist('skip','var')
    skip = 1; % use every nth sequence.
end

% Load Sequence and Counts data. These Sequences are already reverse
% complemented from the reads. So they are 5'V -> 3'J.
data_path = 'data/EXTRACTED_DATA/';
[data, header] = load_data([data_path data_filename]);

%Parse header values
[Cs, Vs] = strtok(header.comments,'=');
Vs=strrep(Vs,'=','');

all_Counts = data.Counts;
N = size(all_Counts,1)
N_read = ceil((N-start+1)/skip)
Read_Length = str2num(Vs{strcmp(Cs,'Read_Length')})
all_Sequences = data.Sequence(start:skip:N);
assert(size(all_Sequences,1)==N_read);
assert(Read_Length==length(all_Sequences{1}));
clear data;

% load genomic segments, as array of structs.
[genV, genD, genJ] = load_Mus_TCR_genomic_segments();

%% Generate output file names.

% Number of sequences processed for each output alignments file.
N_per_file = 5000;

if start==1 && skip==1
    if strcmp(data_filename(end-3:end),'data')
        al_string = [data_filename(1:end-4) 'alignments'];
    end

    % If N is large, split output file into chunks of 150,000.
    n_out_files = ceil(N/N_per_file);
    for i=1:n_out_files
        out_filenames{i} = ['results/' al_string '.' num2str(i) '.mat' ];
    end
else
    out_filenames{1} = ['results/' data_filename '.start.' num2str(start) '.skip.' num2str(skip) '.mat' ];
    n_out_files=1;
end


%% Settings and Initializations.

J_allowed_errors = 1; % Maximum number of errors allowed in J match
V_allowed_errors = 10; % Maximum number of errors allowed in V match

% V and J matches are sorted by score.
% score = match_length - error_cost*n_errors;
% I'm being lenient with errors at this alignment step.
J_error_cost = 4;
V_error_cost = 4;

% Alignment Score Matrix. USED ONLY FOR D MATCHING
scorematrix=nuc44;
scorematrix=scorematrix(1:4,1:4);
scorematrix(scorematrix==-4)=-14;


max_D_aligns = 200;
% Maximum number of deletions from genomic sequence.
J_maximum_deletion = 22;
V_maximum_deletion = 20;
D_maximum_deletion = 16;

% Minimum match length for alignment step. This value can be different in
% the assignments step.
J_minimum_alignment_length = 2;
if Read_Length==101
    V_minimum_alignment_length = 30; % For 101 bp reads, this is 25, for 60 bp reads, this is 15.
else
    V_minimum_alignment_length = 15; % For 101 bp reads, this is 25, for 60 bp reads, this is 15.
end

% To store error and palindrome info for the cases where number of
% deletions is less than the 'min_deletions' (from alignments), how far do
% you go:
negative_excess_deletions_max = 3;


% For saving all of the above parameters to file, I make a struct.
% parfor complains if I refer to the structure fields, so I have above
% variables.
alignment_parameters.scorematrix = scorematrix;
alignment_parameters.J_allowed_errors = J_allowed_errors;
alignment_parameters.V_allowed_errors = V_allowed_errors;
alignment_parameters.J_error_cost = J_error_cost; % used in aligning V and J. score=match_length - error_cost*n_errors;
alignment_parameters.V_error_cost = V_error_cost;
alignment_parameters.J_maximum_deletion = J_maximum_deletion;
alignment_parameters.V_maximum_deletion = V_maximum_deletion;
alignment_parameters.D_maximum_deletion = D_maximum_deletion;
alignment_parameters.J_minimum_alignment_length = J_minimum_alignment_length;
alignment_parameters.V_minimum_alignment_length = V_minimum_alignment_length;
alignment_parameters.negative_excess_deletions_max = negative_excess_deletions_max;
alignment_parameters.read_length = Read_Length;


%% Main loops where alignments are done.

% Loop over output files.
for f=1:n_out_files

    % Take chunk of sequences to be processed for this output file
    N_this_file = min(N_per_file, N_read - (f-1)*N_per_file);
    V = cell(N_this_file,1);
    D = cell(N_this_file,1);
    J = cell(N_this_file,1);

    inds_ok = zeros(N_this_file,1);
    Sequence = all_Sequences(1 + (f-1)*N_per_file : min(f*N_per_file, N_read) );
    Counts = all_Counts(1 + (f-1)*N_per_file : min(f*N_per_file, N_read) );

    % Loop over sequences
    parfor i=1:N_this_file
        if mod(i,100)==0, disp([num2str(i) '***************']); end;
        seq=Sequence{i}; % The sequence is already the rev complement of the raw read.

        % Align with each genomic V. Get list of plausible alignments.
        [seq_v_ok, V{i}] = match_Vs(seq, genV, V_minimum_alignment_length, V_maximum_deletion, negative_excess_deletions_max, V_allowed_errors,  V_error_cost);
        if~seq_v_ok
            %disp(['V align failure for seq ' num2str(i)]);
            continue;
        end

        % Align with each genomic J. Get list of plausible alignments.
        [seq_j_ok, J{i}] = match_Mus_Js(seq, genJ, J_minimum_alignment_length, J_maximum_deletion, negative_excess_deletions_max, J_allowed_errors, J_error_cost);
        if~seq_j_ok
            %disp(['J align failure for seq ' num2str(i)]);
            continue;
        end

        % Align NDN based on longest V and J alignments with each genomic D.
        % Get list of plausible alignments
        V_end =  max(V{i}.align_length); % end of longest V alignment
        J_start = min(J{i}.align_position(1,:)); % start of longest J alignment
        D{i} = match_Ds(seq, genD, V_end, J_start, scorematrix, D_maximum_deletion, negative_excess_deletions_max, max_D_aligns);

        % Store if this sequence was successfully processed.
        inds_ok(i) = 1;
        %disp('Totally fine sequence!')
    end

    % Remove entries that were not successful.
    V(inds_ok==0) = [];
    D(inds_ok==0) = [];
    J(inds_ok==0) = [];
    Sequence(inds_ok==0) = [];
    Counts(inds_ok==0) = [];

    % Convert to regular matrix from cell array. Cell array was used
    % because parfor liked it better.
    V = cell2mat(V);
    D = cell2mat(D);
    J = cell2mat(J);

    ok_fraction = sum(inds_ok)/length(inds_ok);
    % Save this chunk to output file
    save(out_filenames{f},'V','D','J','Sequence','Counts', 'inds_ok','alignment_parameters','ok_fraction');
end

end
