Extract_data_CRISPResso_Clonal_Edited_Cells.m

%Your current folder must contain the main folder 
%containing subfolders of CRISPResso Results.
%Type in command Window the main folder.

prompt = 'Which folder contains subfolders with CRISPResso results? ';
expFolder = input(prompt,'s');
cd(expFolder);
    
%First Matlab moves all html files to a new folder html
    
testhtml = dir('*.html')
testhtml2 = numel(testhtml)

if testhtml2 > 0
mkdir html
movefile *.html html
end

%Matlab defines the list of CRISPResso Results subfolders 
    
FolderList = dir(['.' '/CRISP*']);

%Type if you provided the coding sequence to CRISPResso. 

prompt = 'Did you provide coding sequence to CRISPResso? (yes/no) ';
coding = input(prompt,'s');

%Depending if you provide the coding sequence MATLAB will extract different
%data

switch coding
    
    case 'yes'
  
%In the case of providing the coding sequence.
%The next loop iterates the extraction of the following data
%from each subfolder:

%·Total aligned reads from the CRISPResso_mapping_statistics.txt

%·The inframe and frameshift mutations reads from Frameshift_analysis.txt

%·Total unModified and modified reads from the 
%CRISPResso_quantification_of_editing_frequency.txt

%·Up to tree Allele frequencies and sequences arising from the 
%Alleles_frequency_table_around*.txt

    for i=1:numel(FolderList)

        cd(FolderList(i).name)

    FilenameS = ['CRISPResso_mapping_statistics.txt'];
    FilenameF = ['Frameshift_analysis.txt'];
    FilenameQ = ['CRISPResso_quantification_of_editing_frequency.txt'];

    %Now matlab is going to check if the needed data files are present
    %in the current subfolder. Depending of the CRISPResso analysis succes 
    %there are 2 conditions:

    ExistenceFRAME = exist (FilenameS);
    ExistenceSIMAP = exist (FilenameF);
    TotalExist = ExistenceFRAME + ExistenceSIMAP;

    %1.- Successful analysis when coding sequence is provided 
    %to CRISPResso. All files CRISPResso_mapping_statistics.txt, 
    %Frameshift_analysis.txt and 
    %CRISPResso_quantification_of_editing_frequency.txt are generated by
    %CRISPResso.

    if TotalExist == 4

    %textscan generates variables which
    %contain in their boxes the isolated values of interest.

        fileID = fopen(FilenameF);
        Mutationsarray = textscan(fileID,'%s', 'Delimiter',': ');

        fileID = fopen(FilenameS);
        Readsarray = textscan(fileID,'%s', 'Delimiter','\t');

        fileID = fopen(FilenameQ);
        Editing_Quantification = textscan(fileID,'%s', 'Delimiter','\t');

        variablename = dir ('Alleles_frequency_table_around*.txt');
        Allelesfreqdatatxt = variablename.name;

        Afq = fopen(Allelesfreqdatatxt);
        Allelesfrequencytable = textscan(Afq,'%s', 'Delimiter','\t');

    %Matlab fixes Alignedreads, InFrameMutations, FrameshiftMutations
    %UnModified_reads and UnModified_reads by calling the boxes containing
    %the values of interest.

        Alignedreads{i} = Readsarray{1, 1}{10, 1};
        InFrameMutations{i} = Mutationsarray{1, 1}{9, 1};
        FrameshiftMutations{i} = Mutationsarray{1, 1}{13, 1};
        UnModified_reads{i} = Editing_Quantification{1,1}{26,1};
        Modified_reads{i} = Editing_Quantification{1,1}{27,1};
        Allele1Seq{i} = Allelesfrequencytable{1, 1}{9, 1};
        Allele1Frequency{i} = Allelesfrequencytable{1, 1}{16, 1};
        Ref_Sequence{i} = Allelesfrequencytable{1, 1}{10, 1};

    %Matlab fixes up to tree allele sequences and frequencies in case they
    %are detected by CRISPResso.
                
        switch numel(Allelesfrequencytable{1})
            case 16 %Only one allele is detected
                Allele2Seq{i} = 'NaN';
                Allele2Frequency{i} = 'NaN';
                Allele3Seq{i} = 'NaN';
                Allele3Frequency{i} = 'NaN';
            case 24 %Only two alleles are detected
                Allele2Seq{i} = Allelesfrequencytable{1, 1}{17, 1};
                Allele2Frequency{i} = Allelesfrequencytable{1, 1}{24, 1};
                Allele3Seq{i} = 'NaN';
                Allele3Frequency{i} = 'NaN';                
            otherwise %More than two alleles are detected
                Allele2Seq{i} = Allelesfrequencytable{1, 1}{17, 1};
                Allele2Frequency{i} = Allelesfrequencytable{1, 1}{24, 1}; 
                Allele3Seq{i} = Allelesfrequencytable{1, 1}{25, 1};
                Allele3Frequency{i} = Allelesfrequencytable{1, 1}{32, 1}
        end

    end

    %2.-NOT Successful analysis. Any of the needed files
    %are generated by CRISPResso.

    %Matlab fixes Alignedreads, UnModified_reads, Modified_reads,
    %InFrameMutations and FrameshiftMutations with NaN because in
    %these case are NOT calculated by CRISPResso.

    if TotalExist == 0

        Alignedreads{i} = 'NaN';
        InFrameMutations{i} = 'NaN';
        FrameshiftMutations{i} = 'NaN';
        UnModified_reads{i} = 'NaN';
        Modified_reads{i} = 'NaN';
        Allele1Seq{i} = 'NaN';
        Allele1Frequency{i} = 'NaN';
        Ref_Sequence{i} = 'NaN';
        Allele2Seq{i} = 'NaN';
        Allele2Frequency{i} = 'NaN';
        Allele3Seq{i} = 'NaN';
        Allele3Frequency{i} = 'NaN';

    end


    %Matlab exits the current subfolder to start in the next one.

    fclose('all')

    cd ..

    end

    %Finally Matlab calculates:
    %Percentage of reads: Un-modified, modified, with Frameshift mutations,
    %with In-Frame mutations and with total NHEJ events.
    %Commands as str2num char are used to make variables compatible for 
    %math operations.
    %char transforms cell to character string and str2num 
    %transform to number.

    for i = 1:numel(FolderList)

        FS_percentage(i) = (str2num(char(FrameshiftMutations(i)))...
            /str2num(char(Alignedreads(i)))) * 100;

        IF_percentage(i) = (str2num(char(InFrameMutations(i)))...
            /str2num(char(Alignedreads(i)))) * 100;

        UnMod_percentage(i) = (str2num(char(UnModified_reads(i)))...
            /str2num(char(Alignedreads(i)))) * 100;
        
        Mod_percentage(i) = (str2num(char(Modified_reads(i)))...
            /str2num(char(Alignedreads(i)))) * 100;

    end
    
%Set percentage values cut off for homozigosity or heterozigosity.
%Matlab uses this values to assign genotypes.
%Those samples presenting allele 1 frequency higher than 80% are considered
%homozigotes. Those samples presenting allele 1 and 2 frecuencies
%higher than 40% are considered heterozigotes.

cutoff_hom = 80
cutoff_het = 40
    
    
for i = 1:numel(FolderList)    
    
%wtwt assigment

if (UnMod_percentage(1,i) > cutoff_hom) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_hom)
    wtwt(i) = 1
else
    wtwt(i) = 0
end

%ifwt assigment

if (IF_percentage(1,i) > cutoff_het) & ...
        (UnMod_percentage(1,i) > cutoff_het) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)
    ifwt(i) = 1
else
    ifwt(i) = 0
end

%ifif assigment

if (IF_percentage(1,i) > cutoff_hom) & ...
        (   (str2double(Allele1Frequency(1,i)) > cutoff_hom)  |  ...
        ((str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)))
    ifif(i) = 1
else
    ifif(i) = 0
end
 
%fsif assigment

if (FS_percentage(1,i) > cutoff_het) & ...
        (IF_percentage(1,i) > cutoff_het) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)
    fsif(i) = 1
else
    fsif(i) = 0
end

%fswt assigment

if (FS_percentage(1,i) > cutoff_het) & ...
        (UnMod_percentage(1,i) > cutoff_het) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)
    fswt(i) = 1
else
    fswt(i) = 0
end

%fsfs assigment

if (FS_percentage(1,i) > cutoff_hom) & ...
        (   (str2double(Allele1Frequency(1,i)) > cutoff_hom)  |  ...
        ((str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)))
    fsfs(i) = 1
else
    fsfs(i) = 0
end

%mixed clone assigment

if (wtwt(i) + ifwt(i) + ifif(i) + fsif(i) + fswt(i) + fsfs(i)) == 0
    mixed(i) = 1
else
    mixed(i) = 0
end
  
end
    
%Finally, matlab builds a table that writes everything in a excel file
%(.xlsx file)

    TableMyseqExperiments = table(char(FolderList(1:numel(FolderList)).name),...
        char(Alignedreads{1,:}),...
        char(FrameshiftMutations),...
        char(InFrameMutations), ...
        transpose(UnModified_reads),...
        transpose(FS_percentage), ...
        transpose(IF_percentage), ...
        transpose(UnMod_percentage),...
        transpose(Mod_percentage),...
        char(Allele1Frequency{1,:}),...
        char(Allele1Seq{1,:}),...
        char(Allele2Frequency{1,:}),...
        char(Allele2Seq{1,:}),...
        char(Allele3Frequency{1,:}),...
        char(Allele3Seq{1,:}),...
        char(Ref_Sequence{1,:}), transpose(wtwt), transpose(ifwt),...
        transpose(ifif), transpose(fsif), transpose(fswt),...
        transpose(fsfs), transpose(mixed) )

    TableMyseqExperiments.Properties.VariableNames = {'Sample' ...
        'Aligned_Reads' 'Frameshift_Mutations' 'InFrame_Mutations' ...
        'UnModified_Reads' 'Percentage_FSmut' 'Percentage_IFmut'...
        'Percentage_UnMod' 'Percentage_Mod' 'Allele1_Freq' 'Allele1_Seq' ...
        'Allele2_Freq' 'Allele2_Seq' 'Allele3_Freq' 'Allele3_Seq'...
        'Ref_Sequence' 'wtwt'	'ifwt'	'ifif'	'fsif'	'fswt'	'fsfs'...
        'mixedclone'}


   writetable(TableMyseqExperiments,[expFolder '.xlsx'])
   
    case 'no'
        
%In the case of NOT providing the coding sequence.
%The next loop iterates the extraction of the following data
%from each subfolder:

%·Total aligned reads from the CRISPResso_mapping_statistics.txt

%·Total unModified and Modified reads from the 
%CRISPResso_quantification_of_editing_frequency.txt

%·Up to tree Allele frequencies and sequences arising from the 
%Alleles_frequency_table_around*.txt

    for i=1:numel(FolderList)

        cd(FolderList(i).name)

    FilenameS = ['CRISPResso_mapping_statistics.txt'];
    FilenameQ = ['CRISPResso_quantification_of_editing_frequency.txt'];


    %Now matlab is going to check if the needed data files are present
    %in the current subfolder. Depending of the CRISPResso analysis succes 
    %there are 2 conditions:

    ExistenceFRAME = exist (FilenameS);

    %1.-Successful analysis when coding sequence is NOT
    %provided to CRISPResso. Only CRISPResso_mapping_statistics.txt and 
    %CRISPResso_quantification_of_editing_frequency.txt are generated by
    %CRISPResso.

    if ExistenceFRAME == 2
        
    %textscan command generates variables which
    %contain in their boxes the isolated values of interest.
        
        fileID = fopen(FilenameS);
        Readsarray = textscan(fileID,'%s', 'Delimiter','\t');

        fileID = fopen(FilenameQ);
        Editing_Quantification = textscan(fileID,'%s', 'Delimiter','\t');

        variablename = dir ('Alleles_frequency_table_around*.txt');
        Allelesfreqdatatxt = variablename.name;

        Afq = fopen(Allelesfreqdatatxt);
        Allelesfrequencytable = textscan(Afq,'%s', 'Delimiter','\t');

    %Matlab fixes Alignedreads, UnModified_reads and Modified reads
    %by calling the boxes containing the values of interest.
    %InFrameMutations and FrameshiftMutations 
    %are asigned with NaN because they are NOT calculated by CRISPResso.

        Alignedreads{i} = Readsarray{1, 1}{10, 1};
        UnModified_reads{i} = Editing_Quantification{1,1}{24,1};
        Modified_reads{i} = Editing_Quantification{1,1}{25,1};
        Allele1Seq{i} = Allelesfrequencytable{1, 1}{9, 1};
        Allele1Frequency{i} = Allelesfrequencytable{1, 1}{16, 1};
        Ref_Sequence{i} = Allelesfrequencytable{1, 1}{10, 1};

        switch numel(Allelesfrequencytable{1})
            case 16
                Allele2Seq{i} = 'NaN';
                Allele2Frequency{i} = 'NaN';
                Allele3Seq{i} = 'NaN';
                Allele3Frequency{i} = 'NaN';
            case 24
                Allele2Seq{i} = Allelesfrequencytable{1, 1}{17, 1};
                Allele2Frequency{i} = Allelesfrequencytable{1, 1}{24, 1};
                Allele3Seq{i} = 'NaN';
                Allele3Frequency{i} = 'NaN';                
            otherwise
                Allele2Seq{i} = Allelesfrequencytable{1, 1}{17, 1};
                Allele2Frequency{i} = Allelesfrequencytable{1, 1}{24, 1}; 
                Allele3Seq{i} = Allelesfrequencytable{1, 1}{25, 1};
                Allele3Frequency{i} = Allelesfrequencytable{1, 1}{32, 1}
        end
    end

    %3.-NOT Successful analysis. Any of the needed files
    %are generated by CRISPResso.

    %Matlab fixes Alignedreads, UnModified_reads, 
    %Modified_reads, InFrameMutations and FrameshiftMutations with NaN because in
    %these case are NOT calculated by CRISPResso.

    if TotalExist == 0

        Alignedreads{i} = 'NaN';
        UnModified_reads{i} = 'NaN';
        Modified_reads{i} = 'NaN';
        Allele1Seq{i} = 'NaN';
        Allele1Frequency{i} = 'NaN';
        Ref_Sequence{i} = 'NaN';
        Allele2Seq{i} = 'NaN';
        Allele2Frequency{i} = 'NaN';
        Allele3Seq{i} = 'NaN';
        Allele3Frequency{i} = 'NaN';

    end


    %Matlab exits the current subfolder to start in the next one.

    fclose('all')

    cd ..

    end

    %Finally Matlab calculates:
    %Percentage of reads: Un-modified, with Frameshift mutations, with 
    %In-Frame mutations and with total NHEJ events.
    %Commands as str2num char are used to make variables compatible for 
    %math operations.
    %char transforms cell to character string and str2num transform to number.

    for i = 1:numel(FolderList)

        UnMod_percentage(i) = (str2num(char(UnModified_reads(i)))...
            /str2num(char(Alignedreads(i)))) * 100;
        
        Mod_percentage(i) = (str2num(char(Modified_reads(i)))...
            /str2num(char(Alignedreads(i)))) * 100;

    end
    
%Set percentage values cutoff for homozigosity or heterozigosity

cutoff_hom = 80
cutoff_het = 40
    
    
for i = 1:numel(FolderList)    
    
%wtwt conditional

if (UnMod_percentage(1,i) > cutoff_hom) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_hom)
    wtwt(i) = 1
else
    wtwt(i) = 0
end

%wtmut conditional

if (UnMod_percentage(1,i) > cutoff_het) & ...
        (Mod_percentage(1,i) > cutoff_het) & ...
        (str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)
    wtmut(i) = 1
else
    wtmut(i) = 0
end

%mutmut conditional

if (Mod_percentage(1,i) > cutoff_hom) & ...
        (   (str2double(Allele1Frequency(1,i)) > cutoff_hom)  |  ...
        ((str2double(Allele1Frequency(1,i)) > cutoff_het) & ...
        (str2double(Allele2Frequency(1,i)) > cutoff_het)))
    mutmut(i) = 1
else
    mutmut(i) = 0
end

%mixed clone conditional

if (wtwt(i) + wtmut(i) + mutmut(i)) == 0
    mixed(i) = 1
else
    mixed(i) = 0
end
  
end
    
    
    

    %Finally, matlab builds a table that writes everything in a excel file
    %(.xlsx file)

    TableMyseqExperiments = table(char(FolderList(1:numel(FolderList)).name),...
        char(Alignedreads{1,:}),...
        transpose(UnModified_reads),...
        transpose(Modified_reads),...
        transpose(UnMod_percentage),...
        transpose(Mod_percentage),...
        char(Allele1Frequency{1,:}),...
        char(Allele1Seq{1,:}),...
        char(Allele2Frequency{1,:}),...
        char(Allele2Seq{1,:}),...
        char(Allele3Frequency{1,:}),...
        char(Allele3Seq{1,:}),...
        char(Ref_Sequence{1,:}), transpose(wtwt), transpose(wtmut),...
        transpose(mutmut), transpose(mixed))

    TableMyseqExperiments.Properties.VariableNames = {'Sample' ...
        'Aligned_Reads'  ...
        'UnModified_Reads' 'Modified_Reads'...
        'Percentage_UnMod' 'Percentage_Mod' 'Allele1_Freq' 'Allele1_Seq' ...
        'Allele2_Freq' 'Allele2_Seq' 'Allele3_Freq' 'Allele3_Seq'...
        'Ref_Sequence' 'wtwt'	'wtmut'	'mutmut' ...
        'mixedclone'}


   writetable(TableMyseqExperiments,[expFolder '.xlsx'])

   cd ..
   
end