CRISPResso2parser_clonal-HDR.py

# -*- coding: utf-8 -*-

from __future__ import division
import imp
import glob
import os
import shutil
import pandas as pd
import numpy as np



#Type in command Window the main folder containning your 
#CRISPResso_on.. folders

print("Enter your results folder: ")

expFolder = raw_input()

#Type in command Window the
#Minimun percentage to assume a pure clone. Allele1 percentage must be
#greater than this value to be classified as homozygotes. Allele 1 and
#Allele2 percentages must be higher than this value/2 to be classified as
#heterozygotes. Otherwise clones will be classified as non-pure genotyped clones.


print("Enter your desired clone purity value: ")

Purity_percentage = int(raw_input())

#Type in the command window if you provided coding sequence to CRISPResso2.

print("Have you provided coding sequence to CRISPResso? (YES or NO): ")

CODING_seq = str.upper(raw_input())

#Type in the command window the name of your HDR used
#in your genome editing experiment.

print("Enter your HDR templates names separated by 1 white space: ")

HDR_names = raw_input().split()

#Type in the command window the sequence of your HDR used 
#in your genome editing experiment.

print("Enter your respectively HDR sequences separated by 1 white space: ")
      
HDR_seqs = str.upper(raw_input()).split()

os.chdir(expFolder)


#Moves all html files to a new folder


dirName = 'html'
 
if not os.path.exists(dirName):
    os.mkdir(dirName)
    print("Directory " , dirName ,  " Created ")
else:    
    print("Directory " , dirName ,  " already exists")


if len(glob.glob('*.html'))>0:
    for file in glob.glob('*.html'):
        shutil.move(file,'html')


# glob.glob defines the list of CRISPResso Results subfolders
   
FolderList = glob.glob('CRISP*')

#The next loop iterates the extraction of the following data
#from results files txt that are present in 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

FilenameS = 'CRISPResso_mapping_statistics.txt'
FilenameF = 'Frameshift_analysis.txt'
FilenameQ = 'CRISPResso_quantification_of_editing_frequency.txt'


Aligned_reads = []
InFrameMutations = []
FrameshiftMutations = []
UnModified_reads = []
Modified_reads = []
Allele1_Seq = []
Allele1_Percentage = []
Al1_Ref_Sequence = []
Al1_Indel_value = []
Al1_Unedited = []
Allele2_Seq = []
Allele2_Percentage = []
Al2_Ref_Sequence = []
Al2_Indel_value = []
Al2_Unedited = []
Allele3_Seq = []
Allele3_Percentage = []
Al3_Ref_Sequence = []

for i in range(0,(len(FolderList))):
    
    os.chdir(FolderList[i])
    
    FilenameA = glob.glob('Alleles_frequency_table_around*.txt')
    
#Depending of the CRISPResso analysis output the presence of the 
#results files txt will change.

    # Case1 Successful analysis when coding sequence is provided 
    # to CRISPResso2. All txt files are generated by CRISPResso.

    if (os.path.exists(FilenameF) and os.path.exists(FilenameS) \
        and os.path.exists(FilenameQ)):

        #pd.read_csv reads the txt files and creates variables containing 
        #our values of interest.

        MutationsData = pd.read_csv(FilenameF, sep= ":| ", engine='python')
        
        ReadsData = pd.read_csv(FilenameS, sep="\t")
        
        Editing_Quantification = pd.read_csv(FilenameQ, sep="\t")
        
        Alleles_table = pd.read_csv(FilenameA[0], sep="\t")

        #Append the corresponding values obtained through the for loop.
        
        Aligned_reads.append(ReadsData.iloc[0][2])
        InFrameMutations.append(MutationsData.iloc[1][1])
        FrameshiftMutations.append(MutationsData.iloc[2][1])
        UnModified_reads.append(Editing_Quantification.Unmodified[0])
        Modified_reads.append(Editing_Quantification.Modified[0])
        Allele1_Seq.append(Alleles_table.iloc[0][0])
        Allele1_Percentage.append(Alleles_table.iloc[0][7])
        Al1_Ref_Sequence.append(Alleles_table.iloc[0][1])
        Al1_Indel_value.append(Alleles_table.iloc[0][3]+Alleles_table.iloc[0][4])
        Al1_Unedited.append(str(Alleles_table.iloc[0][2]))
        
        # Fix up to tree allele sequences and percentages in case they
        # are detected by CRISPResso2.  
    
        if (len(Alleles_table) >2): #More than two alleles are detected
        
                Allele2_Seq.append(Alleles_table.iloc[1][0])
                Allele2_Percentage.append(Alleles_table.iloc[1][7])
                Al2_Ref_Sequence.append(Alleles_table.iloc[1][1])
                Al2_Indel_value.append(Alleles_table.iloc[1][3]+Alleles_table.iloc[1][4])
                Al2_Unedited.append(str(Alleles_table.iloc[1][2]))
                Allele3_Seq.append(Alleles_table.iloc[2][0])
                Allele3_Percentage.append(Alleles_table.iloc[2][7])
                Al3_Ref_Sequence.append(Alleles_table.iloc[2][1])
                
        elif (len(Alleles_table) ==1): #Only one allele is detected
        
                Allele2_Seq.append('NaN')
                Allele2_Percentage.append(np.nan)
                Al2_Ref_Sequence.append('NaN')
                Al2_Indel_value.append(np.nan)
                Al2_Unedited.append('NaN')
                Allele3_Seq.append('NaN')
                Allele3_Percentage.append(np.nan)
                Al3_Ref_Sequence.append('NaN')
        
        elif (len(Alleles_table) ==2):# Only two alleles are detected
        
                Allele2_Seq.append(Alleles_table.iloc[1][0])
                Allele2_Percentage.append(Alleles_table.iloc[1][7])
                Al2_Ref_Sequence.append(Alleles_table.iloc[1][1])
                Al2_Indel_value.append(Alleles_table.iloc[1][3]+Alleles_table.iloc[1][4])
                Al2_Unedited.append(str(Alleles_table.iloc[1][2]))
                Allele3_Seq.append('NaN')
                Allele3_Percentage.append(np.nan)
                Al3_Ref_Sequence.append('NaN')
                
        os.chdir("..")
    
    #Case2 Successful analysis when coding sequence is NOT provided to 
    #CRISPResso2. Only CRISPResso_mapping_statistics.txt and 
    #CRISPResso_quantification_of_editing_frequency.txt are generated by
    #CRISPResso2.

    elif (not os.path.exists(FilenameF) and os.path.exists(FilenameS) \
        and os.path.exists(FilenameQ)):

        ReadsData = pd.read_csv(FilenameS, sep="\t")
        
        Editing_Quantification = pd.read_csv(FilenameQ, sep="\t")

    #Append the corresponding values obtained through the for loop.
    #InFrameMutations and FrameshiftMutations are appended with NaN because 
    #they are NOT calculated by CRISPResso2.

        Aligned_reads.append(ReadsData.iloc[0][2])
        InFrameMutations.append(np.nan)
        FrameshiftMutations.append(np.nan)
        UnModified_reads.append(Editing_Quantification.Unmodified[0])
        Modified_reads.append(Editing_Quantification.Modified[0])
        Allele1_Seq.append(Alleles_table.iloc[0][0])
        Allele1_Percentage.append(Alleles_table.iloc[0][7])
        Al1_Ref_Sequence.append(Alleles_table.iloc[0][1])
        Al1_Indel_value.append(Alleles_table.iloc[0][3] + Alleles_table.iloc[0][4])
        Al1_Unedited.append(str(Alleles_table.iloc[0][2]))
        
        # Fixe up to tree allele sequences and percentages in case they
        # are detected by CRISPResso2.    
    
        if (len(Alleles_table) >2): #More than two alleles are detected
        
                Allele2_Seq.append(Alleles_table.iloc[1][0])
                Allele2_Percentage.append(Alleles_table.iloc[1][7])
                Al2_Ref_Sequence.append(Alleles_table.iloc[1][1])
                Al2_Indel_value.append(Alleles_table.iloc[1][3] + Alleles_table.iloc[1][4])
                Al2_Unedited.append(str(Alleles_table.iloc[1][2]))
                Allele3_Seq.append(Alleles_table.iloc[2][0])
                Allele3_Percentage.append(Alleles_table.iloc[2][7])
                Al3_Ref_Sequence.append(Alleles_table.iloc[2][1])
                
        elif (len(Alleles_table) ==1): #Only one allele is detected
        
                Allele2_Seq.append('NaN')
                Allele2_Percentage.append(np.nan)
                Al2_Ref_Sequence.append('NaN')
                Al2_Indel_value.append(np.nan)
                Al2_Unedited.append('NaN')
                Allele3_Seq.append('NaN')
                Allele3_Percentage.append(np.nan)
                Al3_Ref_Sequence.append('NaN')
        
        elif (len(Alleles_table) ==2):# Only two alleles are detected
        
                Allele2_Seq.append(Alleles_table.iloc[1][0])
                Allele2_Percentage.append(Alleles_table.iloc[1][7])
                Al2_Ref_Sequence.append(Alleles_table.iloc[1][1])
                Al2_Indel_value.append(Alleles_table.iloc[1][3]+Alleles_table.iloc[1][4])
                Al2_Unedited.append(str(Alleles_table.iloc[1][2]))
                Allele3_Seq.append('NaN')
                Allele3_Percentage.append(np.nan)
                Al3_Ref_Sequence.append('NaN')
                
        os.chdir("..")

    #Case3 NOT Successful analysis. Any of the needed files are generated by 
    #CRISPResso2. Python appends Alignedreads, UnModified_reads, 
    #Modified_reads,InFrameMutations and FrameshiftMutations with NaN because
    #in these case are NOT calculated by CRISPResso.

    else:
    
        Aligned_reads.append(np.nan)
        InFrameMutations.append(np.nan)
        FrameshiftMutations.append(np.nan)
        UnModified_reads.append(np.nan)
        Modified_reads.append(np.nan)
        Allele1_Seq.append('NaN')
        Allele1_Percentage.append(np.nan)
        Al1_Ref_Sequence.append('NaN')
        Al1_Indel_value.append(np.nan)
        Al1_Unedited.append('NaN')
        Allele2_Seq.append('NaN')
        Allele2_Percentage.append(np.nan)
        Al2_Ref_Sequence.append('NaN')
        Al2_Indel_value.append(np.nan)
        Al2_Unedited.append('NaN')
        Allele3_Seq.append('NaN')
        Allele3_Percentage.append(np.nan)
        Al3_Ref_Sequence.append('NaN')
    
        os.chdir("..")

#Calculate Percentage of reads: Un-modified, modified, 
#with Frameshift mutations and with In-Frame mutations.

FS_percentage = []
IF_percentage = []
UnMod_percentage = []
Mod_percentage = []


for i in range(0,(len(FolderList))):
    
    FS_percentage.append\
        (FrameshiftMutations[i] / Aligned_reads[i]*100)
    IF_percentage.append\
        (InFrameMutations[i] / Aligned_reads[i]*100)
    UnMod_percentage.append\
        (UnModified_reads[i] / Aligned_reads[i]*100)
    Mod_percentage.append\
        (Modified_reads[i] / Aligned_reads[i]*100)
        
#Generate outputs of Alleles to REF Aligments:

    Al1andREF = zip(Al1_Ref_Sequence, Allele1_Seq)
    Aligment_Allele1 = ['\n'.join(temp) for temp in Al1andREF]
    
    Al2andREF = zip(Al2_Ref_Sequence, Allele2_Seq)
    Aligment_Allele2 = ['\n'.join(temp) for temp in Al2andREF]
    
    Al3andREF = zip(Al3_Ref_Sequence, Allele3_Seq)
    Aligment_Allele3 = ['\n'.join(temp) for temp in Al3andREF]

#Determine Zygosity

Zygosity = []

for i in range(0,(len(FolderList))):
    
    #Homozygote assigment
    
    if (Allele1_Percentage[i] > Purity_percentage):
        
        Zygosity.append('HOMOZYGOTE')
        
    #Heterozygote assigment
        
    elif (Allele1_Percentage[i] > Purity_percentage / 2 and \
                Allele2_Percentage[i] > Purity_percentage / 2):
        
        Zygosity.append('HETEROZYGOTE')
    
    else: 
        
        Zygosity.append('Non-pure')


Genotype_Allele1 = range(len(FolderList))
Genotype_Allele2 = range(len(FolderList))

#Determine Genotype of Allele1 and Allele2 when Coding Seguence is provided

if (CODING_seq == 'YES'):
    
    #Determine Genotype Allele1
    
    for i in range(0,(len(Genotype_Allele1))): 
        
            for j in range(0, (len(HDR_names))):
        
            #Non-pure
            
                if (Zygosity[i] == 'Non-pure'):
                
                    Genotype_Allele1[i] = 'Non-assigned'
                
            #HDR assigment
        
                
                elif (Allele1_Percentage[i] > Purity_percentage / 2 and \
                    Allele1_Seq[i] in HDR_seqs[j]):
                    
                    Genotype_Allele1[i] = HDR_names[j]
            
            #WT assigment
            
                elif (Genotype_Allele1[i] == i and \
                      Al1_Unedited[i] == 'True'):
                    
                    Genotype_Allele1[i] = 'wt'
                    
            #if assigment
            
                elif (Genotype_Allele1[i] == i and \
                    Al1_Unedited[i] == 'False' and \
                    Al1_Indel_value[i] % 3 == 0):
                    
                    Genotype_Allele1[i] = 'if'
                    
            #FS assigment
            
                elif (Genotype_Allele1[i] == i and \
                    Al1_Unedited[i] == 'False' and \
                    Al1_Indel_value[i] % 3 <> 0):

                    Genotype_Allele1[i] = 'fs'
                    
    #Determine Genotype Allele2
    
    for i in range(0,(len(FolderList))): 
        
    
            for j in range(0, (len(HDR_names))):
        
            #Non-assigned because homozygote or non-pure
            
                if (Zygosity[i] == 'Non-pure' or Zygosity[i] == 'HOMOZYGOTE'):
                
                    Genotype_Allele2[i] = 'Non-assigned'
                
            #HDR assigment
                
                elif (Allele2_Percentage[i] > Purity_percentage / 2 and \
                    Allele2_Seq[i] in HDR_seqs[j]):
                    
                    Genotype_Allele2[i] = HDR_names[j]
            
            #WT assigment
            
                elif (Genotype_Allele2[i] == i and \
                    Al2_Unedited[i] == 'True'):
                    
                    Genotype_Allele2[i] = 'wt'
                    
            #if assigment
            
                elif (Genotype_Allele2[i] == i and \
                    Al2_Unedited[i] == 'False' and \
                    Al2_Indel_value[i] % 3 == 0):
                    
                    Genotype_Allele2[i] = 'if'
                    
            #FS assigment
            
                elif (Genotype_Allele2[i] == i and \
                    Al2_Unedited[i] == 'False' and \
                    Al2_Indel_value[i] % 3 <> 0):

                    Genotype_Allele2[i] = 'fs'

#Determine Genotype Allele1 and Allele2 when Coding Seguence is NOT provided

if (CODING_seq == 'NO'):
    
    #Determine Genotype Allele1
    
    for i in range(0,(len(FolderList))): 
        
            for j in range(0, (len(HDR_names))):
        
            #Non-pure
            
                if (Zygosity[i] == 'Non-pure'):
                
                    Genotype_Allele1[i] = 'Non-assigned'
                
            #HDR assigment
        
                
                elif (Allele1_Percentage[i] > Purity_percentage / 2 and \
                    Allele1_Seq[i] in HDR_seqs[j]):
                    
                    Genotype_Allele1[i] = HDR_names[j]
            
            #WT assigment
            
                elif (Genotype_Allele1[i] == i and \
                    Al1_Unedited[i] == 'True'):
                    
                    Genotype_Allele1[i] = 'wt'
                    
            #Mut assigment
            
                elif (Genotype_Allele1[i] == i and \
                    Al1_Unedited[i] == 'False'):
                    
                    Genotype_Allele1[i] = 'mut'

                    
    #Determine Genotype Allele2
    
    for i in range(0,(len(FolderList))): 
        
    
            for j in range(0, (len(HDR_names))):
        
            #Non-assigned because homozygote or non-pure
            
                if (Zygosity[i] == 'Non-pure' or Zygosity[i] == 'HOMOZYGOTE'):
                
                    Genotype_Allele2[i] = 'Non-assigned'
                
            #HDR assigment
                
                elif (Allele2_Percentage[i] > Purity_percentage / 2 and \
                    Allele2_Seq[i] in HDR_seqs[j]):
                    
                    Genotype_Allele2[i] = HDR_names[j]
            
            #WT assigment
            
                elif (Genotype_Allele2[i] == i and \
                    Al2_Unedited[i] == 'True'):
                    
                    Genotype_Allele2[i] = 'wt'
                    
            #mut assigment
            
                elif (Genotype_Allele2[i] == i and \
                    Al2_Unedited[i] == 'False'):
                    
                    Genotype_Allele2[i] = 'mut'
                    
#Build the excel file if Coding sequence is provided.

if (CODING_seq == 'YES'):

    df = pd.DataFrame({'Folder_Sample' : FolderList, \
    'Aligned_reads' : Aligned_reads, \
    'UnModified_reads' : UnModified_reads, \
    'Modified_reads' : Modified_reads, \
    'InFrameMutations' : InFrameMutations, \
    'FrameshiftMutations' : FrameshiftMutations, \
    'Mod_percentage' : Mod_percentage, \
    'UnMod_percentage' : UnMod_percentage, \
    'FS_percentage' : FS_percentage, \
    'IF_percentage' : IF_percentage, \
    'Allele1_Percentage' : Allele1_Percentage, \
    'Aligment_REF_Allele1' : Aligment_Allele1, \
    'Allele2_Percentage' : Allele2_Percentage, \
    'Aligment_REF_Allele2' : Aligment_Allele2, \
    'Allele3_Percentage' : Allele3_Percentage, \
    'Aligment_REF_Allele3' : Aligment_Allele3, \
    'Zygosity' : Zygosity, \
    'Genotype_Allele1' : Genotype_Allele1,\
    'Genotype_Allele2' : Genotype_Allele2}) 
        
    #Establish the desired order of columns
    
    df = df[['Folder_Sample',\
     'Aligned_reads',\
     'Modified_reads',\
     'UnModified_reads',\
     'InFrameMutations',\
     'FrameshiftMutations',\
     'Mod_percentage',\
     'UnMod_percentage',\
     'FS_percentage',\
     'IF_percentage',\
     'Allele1_Percentage', \
     'Aligment_REF_Allele1', \
     'Allele2_Percentage', \
     'Aligment_REF_Allele2', \
     'Allele3_Percentage', \
     'Aligment_REF_Allele3', \
     'Zygosity', \
     'Genotype_Allele1', \
     'Genotype_Allele2' ]]
        
    os.chdir("..")    
    
    #Export to xlsx file
    
    try:
        imp.find_module('openpyxl')
        df.to_excel\
            (expFolder + '_HDR-wCDS.xlsx', float_format="%.2f",na_rep='NaN', index = None)
    except ImportError:
            print('If you want an output excel file Install the module openpyxl'\
                  ' module by typing\n"pip install openpyxl" '\
                      'in Terminal and run again the script.')

#Build the excel file if Coding sequence is NOT provided.

if (CODING_seq == 'NO'):

    df = pd.DataFrame({'Folder_Sample' : FolderList, \
    'Aligned_reads' : Aligned_reads, \
    'UnModified_reads' : UnModified_reads, \
    'Modified_reads' : Modified_reads, \
    'Mod_percentage' : Mod_percentage, \
    'UnMod_percentage' : UnMod_percentage, \
    'Allele1_Percentage' : Allele1_Percentage, \
    'Aligment_REF_Allele1' : Aligment_Allele1, \
    'Allele2_Percentage' : Allele2_Percentage, \
    'Aligment_REF_Allele2' : Aligment_Allele2, \
    'Allele3_Percentage' : Allele3_Percentage, \
    'Aligment_REF_Allele3' : Aligment_Allele3, \
    'Zygosity' : Zygosity, \
    'Genotype_Allele1' : Genotype_Allele1,\
    'Genotype_Allele2' : Genotype_Allele2}) 
        
    #Establish the desired order of columns
    
    df = df[['Folder_Sample',\
     'Aligned_reads',\
     'Modified_reads',\
     'UnModified_reads',\
     'Mod_percentage',\
     'UnMod_percentage',\
     'Allele1_Percentage', \
     'Aligment_REF_Allele1', \
     'Allele2_Percentage', \
     'Aligment_REF_Allele2', \
     'Allele3_Percentage', \
     'Aligment_REF_Allele3', \
     'Zygosity', \
     'Genotype_Allele1', \
     'Genotype_Allele2' ]]
        
    os.chdir("..")    
    
    #Export to xlsx file
    
    try:
        imp.find_module('openpyxl')
        df.to_excel\
            (expFolder + '_HDR-woCDS.xlsx', float_format="%.2f",na_rep='NaN', index = None)
    except ImportError:
            print('If you want an output excel file Install the module openpyxl'\
                  ' module by typing\n"pip install openpyxl" '\
                      'in Terminal and run again the script.')
            
print('\n\nThank you for using CRISPResso2parser_clonal-HDR.py \n\n'\
      'Enjoy the view of your CRISPResso2 results summary! :)\n\n')