Skip to content
Snippets Groups Projects
Commit 958e2857 authored by Rohmer Coralie's avatar Rohmer Coralie
Browse files

add rate subsitution, insertion, deletion in data

parent 49314ae9
No related branches found
No related tags found
No related merge requests found
......@@ -73,6 +73,8 @@ for file in files_alignement:
occurence_similarity = Counter(similarity)
nb_ambiguity=0
nb_identity=0
nb_deletion=0
nb_insertion=0
for i in "RYSWKMBDHVN":
if (i in occurence_seq):
......@@ -88,6 +90,12 @@ for file in files_alignement:
for i in range(0,len(seq)):
if(seq[i]==hsp.hit[i]):
nb_identity+=1
else:
if(seq[i] == "-" ):
nb_deletion+=1
else:
if(hsp.hit[i]== "-" ):
nb_insertion+=1
if (MSA not in data):
data[MSA]={}
......@@ -95,12 +103,16 @@ for file in files_alignement:
data[MSA][read_size] = {}
match=round(tool_consensus.score_between_two_iupac(hsp.hit.seq,seq.seq),1)
error=round(len(seq)-match,1)
nb_substitution = round(error - nb_deletion - nb_insertion,1)
r_ambiguity=round( (nb_ambiguity/(size_seq)*100) ,1)
r_identity=round( (nb_identity/len(seq)*100) ,1)
r_match=round( (match/len(seq)*100) ,1)
r_error=round( (error/len(seq)*100) ,1)
data[MSA][read_size][nb_read]=[nb_ambiguity,r_ambiguity,nb_identity,r_identity,
error,r_error,match,r_match,size_seq]
r_sub=round((nb_substitution/error*100),1)
r_del=round((nb_deletion/error*100),1)
r_ins=round((nb_insertion/error*100),1)
data[MSA][read_size][nb_read]=[nb_ambiguity,r_ambiguity,nb_identity,r_identity,match,r_match,
error,r_error,nb_substitution,r_sub,nb_deletion,r_del,nb_insertion,r_ins,size_seq]
if (TIME == True):
pass
......@@ -125,7 +137,7 @@ if (TIME == True):
#output
sep=","
print("MSA","region_size","depth","number_Ambiguity","percentage_Ambiguity","number_Identity","percentage_Identity","number_Error","percentage_Error","number_Match","percentage_Match","size",sep=sep,end="")
print("MSA","region_size","depth","number_Ambiguity","percentage_Ambiguity","number_Identity","percentage_Identity","number_Match","percentage_Match","number_Error","percentage_Error","number_Substitution","percentage_Substitution","number_Deletion","percentage_Deletion","number_Insertion","percentage_Insertion","size",sep=sep,end="")
if (TIME == True):
print(sep+"time","elapsed","memory",sep=sep)
else:
......
......@@ -196,7 +196,7 @@ def consensus(alignment,threshold):
seq+="-"
return seq
#Creation of the consensus sequence
#Index jaccard for two iupac
def index_jaccard(nuc_ref,nuc_cons):
nb_nuc_ref = len(nuc_ref)
nb_nuc_cons = len(nuc_cons)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment