#!/usr/bin/python3 # -*- coding: utf-8 -*- ############################################### #Dresscode : class : uppercase # function : begins with a capital # variable : words separated by capital # object : words separated by underscore from functionObjectVCF_creator import * import re import os import sys import subprocess import re import time # # CheckStrandAndReverseNucleotide(self,nucleo):"""Reverse the alt nucleotide if it is needed""" # char2char = dict() # for fast reverse complement computations # char2char['A'] = 'T' # char2char['T'] = 'A' # char2char['C'] = 'G' # char2char['G'] = 'C' # char2char['a'] = 't' # char2char['t'] = 'a' # char2char['c'] = 'g' # char2char['g'] = 'c' # def ReverseComplement(nucleotide): # """Take a sequence or a nucleotide and reverse it""" # return ''.join(char2char[c] for c in nucleotide)[::-1] rev = str.maketrans("acgtACGT", "tgcaTGCA") def ReverseComplement(seq: str) -> str: return seq.translate(rev)[::-1] #Class and methods_________________________________________________________________________________________________________ #INDEX_____________________________________________________________________________________________________________ #________________class VARIANT(): corresponds to a discosnp bubble with information common to both paths________________ # ______ #__/ \__ # \______/ # => Methods # FillInformationFromHeader(self,VCFObject):"""Parsing of the DiscoSnp++ header""" # CheckContigUnitig(self,unitig,contig):"""Checks if there is an extension in the form of contig or unitig to take into account in the mapping position""" # RetrievePolymorphismFromHeader(self):'''Gets from the dicoAllele all the positions, and the nucleotides''' # MismatchChecker(self):"""In case of divergent main position (case snp whose two paths are mapped ) to define the reference = > check the number of mismatch # ( If the number of mismatch is the same in both cases it is the lower lexicographical SNP which is selected for reference . # The Boolean allows to know the reference SNP ) """ # RetrieveGenotypes(self,nbGeno,VCFObject):"""Gets the genotype, the coverage and the likelihood by sample and print it in the correspondand fields. The genotype is determined by DiscoSnp++ (which considered the upper path as reference). If the “REF” corresponds the upper path, the genotype in the VCF is identical to the genotype in DiscoSnp++, else it's the opposite ( 1/1 becomes 0/0 and so on).""" # FillVCF(self,VCFfile,nbGeno,table,VCFObject): """Take all necessary input variables to fill the vcf; Fills the fields of the table which will be printed in the vcf ; return the table""" # WhichPathIsTheRef(self,VCFObject): """Finds which path is identical to the reference genome and defines it as the ref : specific method for each type of variant""" # RetrieveMappingPositionCouple(self): """Defines the mapping position for the couple of variant""" # CheckCoupleVariantID(self):"""Test if the couple of paths has the same ID""" #________________class PATH(): """corresponds to one path"""________________ # RetrieveSeq(self,seq):"""Getter for sequence""" # RetrieveDicoMappingPosition(self):"""Retrieves for each path alignment information in a list ; retrieves a dictionary with all the positions of a path and the number of associated mismatch""" # CheckBitwiseFlag(self):"""Checks if the BitwiseFlag contains the tested value such as : read reverse strand, read unmmaped and so on.""" # CigarcodeChecker(self):"""Checks in the cigarcode of the samfile if there is a shift in the alignment between the path and the reference""" # ReferenceChecker(self,shift,posCentraleRef,VCFObject):"""Function which allows to get the MD tag parsing; checks if path nucleotide is identical to the reference nucleotide""" # RetrieveCoverage(self):"""Get the coverage by path in the discosnp++ header""" # GetTag(self):"""Gets the number of mismatch in the samline""" # CheckPosVariantFromRef(self,VCFObject): """Checks if the variant is identical to the reference or not ; defines the nucleotide on the reference""" #________________class SNP(VARIANT):________________ # WhichPathIsTheRef(self,VCFObject):"""""" #________________class INDEL(VARIANT):________________ # RetrievePolymorphismFromHeader(self):"""""" # WhichPathIsTheRef(self,VCFObject):"""""" # RetrieveMappingPositionCouple(self):"""""" #________________class SNPSCLOSE(VARIANT):________________ # RetrieveDicoClose(self,dicoCloseUp,dicoCloseLow): # WhichPathIsTheRef(self,VCFObject): # FillVCF(self,VCFfile,nbGeno,table,VCFObject): #________________class VCFFIELD():________________ # PrintOneLine(self,table,VCF): def shift_from_cigar_code(cigarcode, pospol): # print ("shift",cigarcode,pospol) parsingCigarCode=re.findall(r'(\d+|[A-Za-z])',cigarcode) #ParsingCigarCode=['2', 'S', '3', 'M', '1', 'I', '25', 'M'] # print (parsingCigarCode) shift=0 pos=0 i=1 while i=pospol: # print("return", str(pospol+shift)) return pospol+shift#takes into account the shift to add the after the mapping position and the variant position in the sequence i+=2 # print("return", str(pospol+shift)) return pospol+shift#takes into account the shift to add the after the mapping position and the variant position in the sequence class VARIANT(): """Object corresponding to a discosnp++ bubble""" def __init__(self,line1,line2): self.upper_path = PATH(line1)#line in the file corresponding to the upper path self.lower_path = PATH(line2)#line in the file corresponding to the lower path self.variantID = "" #ID of discoSnp++ self.discoName = "" #name of the variant : SNP_higher_path_99 self.len_unitig_left = 0#length of the unitig left self.len_unitig_right = 0#length of the unitig right self.len_contig_left = 0#length of the contig left self.len_contig_right = 0#length of the contig right self.rank = ""#rank calculated by discosnp++ self.nb_pol = ""# number of polymorphism in the disco path self.dicoGeno = {}#dictionnary with the information of the genotype dicoGeno[listgeno[0]] = [listgeno[1],listlikelihood] self.dicoAllele = {}#dictionnary of all the information from the header of discosnp++ : depending on the variant self.mappingPositionCouple = 0#mapping position of the bubble after correction self.dicoIndex = {}#dictionnary with all the index of every items in discoSnp++ header #--------------------------------------------------------------------------------------------------------------------------- ##Example :SNP_higher_path_99|P_1:30_A/C|high|nb_pol_1|left_unitig_length_129|right_unitig_length_901|C1_0|C2_30|G1_1/1:744,116,6|G2_0/0:6,95,604|rank_1.00000 #--------------------------------------------------------------------------------------------------------------------------- def FillInformationFromHeader(self,VCFObject): """Parsing of the DiscoSnp++ header. Gets unitig, contig, rank, genotypes""" headerVariantUp = self.upper_path.listSam[0]#header of the upper path # headerVariantLow= self.lower_path.listSam[0]#header of the lower path discoList=headerVariantUp.split('|')#splitting the header of discosnp++ into a list self.discoName=discoList[0]#fills the attribut discoName of the variant object # print discoList self.variantID = self.discoName.split("_")[-1]#Gets the variantID ex:56468 if self.discoName.startswith("SNP"): VCFObject.variantType="SNP" else: VCFObject.variantType="INDEL" # VCFObject.variantType = self.discoName.split("_")[0]#fills the VCF object with the variant type listgeno=[]#splitted informations by genotype contained in the header of discoSnp++ #Get dicoAllele P_1:30_A/G => {'P_1': ['30', 'A', 'G']} pos=discoList[self.dicoIndex["P_"]].split(',') if VCFObject.variantType=="SNP":#Specific dictionary dicoAllele in case of simple snp for item in pos: listP_=item.split(":") self.dicoAllele[listP_[0]]=[(listP_[1].split("_")[0]),(listP_[1].split("_")[1].split("/")[0]),(listP_[1].split("_")[1].split("/")[1])]#{'P_1': ['30', 'A', 'G']} elif VCFObject.variantType=="INDEL":#INDEL case : P_1:30_3_2 => {'P_1': ['30', '3', '2']} specific to the INDEL for item in pos: listP_=item.split(":") self.dicoAllele[listP_[0]]=[(listP_[1].split("_")[0]),(listP_[1].split("_")[1]),(listP_[1].split("_")[2])]# {'P_1': ['30', '3', '2']} # for item in pos: # if VCFObject.variantType=="SNP":#Specific dictionary dicoAllele in case of simple snp # # if ":" in item: # listP_=item.split(":") # self.dicoAllele[listP_[0]]=[(listP_[1].split("_")[0]),(listP_[1].split("_")[1].split("/")[0]),(listP_[1].split("_")[1].split("/")[1])]#{'P_1': ['30', 'A', 'G']} # elif VCFObject.variantType=="INDEL":#INDEL case : P_1:30_3_2 => {'P_1': ['30', '3', '2']} specific to the INDEL # listP_=item.split(":") # self.dicoAllele[listP_[0]]=[(listP_[1].split("_")[0]),(listP_[1].split("_")[1]),(listP_[1].split("_")[2])]# {'P_1': ['30', '3', '2']} #Get unitig if "unitig" in self.dicoIndex and "unitig" in headerVariantUp: self.len_unitig_left=int(discoList[self.dicoIndex["unitig"][0]].split("_")[3]) self.len_unitig_right=int(discoList[self.dicoIndex["unitig"][1]].split("_")[3]) #Get contig if "contig" in self.dicoIndex and "contig" in headerVariantUp: self.len_contig_left=int(discoList[self.dicoIndex["contig"][0]].split("_")[3]) self.len_contig_right=int(discoList[self.dicoIndex["contig"][1]].split("_")[3]) #Get rank if "rank" in self.dicoIndex and "rank" in headerVariantUp: self.rank=discoList[self.dicoIndex["rank"]].split('_')[1] #Get nb_pol self.nb_pol=int(discoList[self.dicoIndex["nb_pol"]].split('_')[2]) if "G" in self.dicoIndex and "G1" in headerVariantUp: for i in self.dicoIndex["G"]: listgeno=discoList[i].replace("_",":").split(":") if len(listgeno)>2: self.dicoGeno[listgeno[0]]=[listgeno[1],(listgeno[2].split(','))] else: self.dicoGeno[listgeno[0]]=[listgeno[1]] #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def RetrievePolymorphismFromHeader(self): '''Gets from the dicoAllele all the positions, and the nucleotides for each variant ''' for key,(posD,ntUp,ntLow) in list(self.dicoAllele.items()): #Goes through the dictionary of parsed header self.upper_path.listPosForward.append(int(posD)+1) self.lower_path.listPosForward.append(int(posD)+1) self.upper_path.listPosReverse.append(len(self.upper_path.seq)-int(posD)) self.lower_path.listPosReverse.append(len(self.upper_path.seq)-int(posD)) self.upper_path.listNucleotideForward.append(ntUp) self.lower_path.listNucleotideForward.append(ntLow) self.upper_path.listNucleotideReverse.append(ReverseComplement(ntUp)) self.lower_path.listNucleotideReverse.append(ReverseComplement(ntLow)) #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def MismatchChecker(self): """In case of divergent main position (case snp whose two paths are mapped ) to define the reference = > check the number of mismatch ( If the number of mismatch is the same in both cases it is the lower lexicographical SNP which is selected for reference . The Boolean allows to know the reference SNP : It fills boolRef ) """ nmUp=None nmLow=None #Two paths mapped if self.upper_path.mappingPosition>0 and self.lower_path.mappingPosition>0: #Checks if both paths are mapped nmUp=int(self.upper_path.dicoMappingPos[self.upper_path.mappingPosition][0]) #Distance with the reference for the snpUp nmLow=int(self.lower_path.dicoMappingPos[self.lower_path.mappingPosition][0]) #Distance with the reference for the snpLow # Pierre 25 06 2021: BUG here in case the two paths are not mapped in the same orientation. if nmUpnmLow : #Checks if the lower path has a distance with the reference smaller than the upper path self.lower_path.boolRef=True self.upper_path.boolRef=False self.upper_path.nucleoRef = self.lower_path.nucleoRef elif nmUp==nmLow: #Checks if both path have the same number of difference if self.discoName.split("_")[0]!="INDEL": #In case of simple snp if self.upper_path.nucleo self.lower_path.nucleo: #Checks the lexicographical order self.lower_path.boolRef=True self.upper_path.boolRef=False self.upper_path.nucleoRef = self.lower_path.nucleoRef else : #If none of the alleles is lexicographically less : checks the mapping position and keeps the lefmost position if self.upper_path.mappingPosition0 and \ self.lower_path.mappingPosition>0) \ and self.upper_path.mappingPosition != self.lower_path.mappingPosition): self.MismatchChecker() return if self.upper_path.boolRef == self.lower_path.boolRef: self.MismatchChecker() return #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def RetrieveMappingPositionCouple(self): #Validation SNP second part (specific method for close snps) """Defines the mapping position for the couple of variant by checking boolRef""" #for INDEL and simple snp if self.upper_path.boolRef==True: self.mappingPositionCouple = self.upper_path.mappingPosition+int(self.upper_path.correctedPos[0])+int(self.mappingPositionCouple)-1 else: self.mappingPositionCouple = self.lower_path.mappingPosition+int(self.lower_path.correctedPos[0])+int(self.mappingPositionCouple)-1 #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- # def CheckStrandAndReverseNucleotideNONOPTIMIZED(self,nucleo): # """Reverse the alt nucleotide if it is needed""" # if self.upper_path.boolRef==True:#Checks if the upper path is the reference # if self.upper_path.boolReverse= = self.lower_path.boolReverse :#if the mapping strand is the same on both path => returns the nucleotide # return(nucleo) # if self.upper_path.boolReverse=="1" and self.lower_path.boolReverse==".": # return (nucleo) # if self.upper_path.boolReverse! = self.lower_path.boolReverse:#if the mapping strand is different on both path => returns the reverse nuclotide # return (ReverseComplement(nucleo)) # elif self.lower_path.boolRef==True:#Checks if the lower path is the reference # if self.upper_path.boolReverse= = self.lower_path.boolReverse or (self.lower_path.boolReverse==1 and self.upper_path.boolReverse=="."):#if the mapping strand is the same on both path => returns the nucleotide # return (nucleo) # if self.lower_path.boolReverse=="1" and self.upper_path.boolReverse==".": # return (nucleo) # if self.upper_path.boolReverse! = self.lower_path.boolReverse:#if the mapping strand is different on both path => returns the reverse nucleotide # return (ReverseComplement(nucleo)) # else : # return (nucleo) # #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def CheckStrandAndReverseNucleotide(self,nucleo): """Reverse the alt nucleotide if it is needed""" if self.upper_path.boolReverse != self.lower_path.boolReverse:# if the mapping strand is different on both path => returns the reverse nuclotide if (str(self.upper_path.boolReverse) == "1" and self.lower_path.boolReverse==".") or (self.upper_path.boolReverse=="." and str(self.lower_path.boolReverse)=="1"): return nucleo # Case in which one of the two variants matched with a 2 or more indel, thus considered as non mapped. June 3017 return ReverseComplement(nucleo) else : return nucleo #--------------------------------------------------------------------------------------------------------------------------- #Example of supplementary alignment #INDEL_higher_path_17964|P_1:30_10_8|low|nb_pol_1|left_unitig_length_346|right_unitig_length_815|C1_12|C2_1|G1_0/1:321,17,162|G2_1/1:848,120,10|rank_0.46189 0 gi|224384768|gb|CM000663.1| 191102952 60 52M * 0 0 AAGAAAAAAGAAATAAAAAAAGAAAAAAAAACGAAATAGCCAGAAGGAATGA * NM:i:2 MD:Z:1G9C40 AS:i:45 XS:i:23 #INDEL_lower_path_17964|P_1:30_10_8|low|nb_pol_1|left_unitig_length_346|right_unitig_length_815|C1_20|C2_43|G1_0/1:321,17,162|G2_1/1:848,120,10|rank_0.46189 0 gi|224384768|gb|CM000663.1| 191102966 123S8M1I30M * 0 0 AAGAAAAAAGAAATAAAAAAAGAAAAAAAAGAAAAAAAAAACGAAATAGCCAGAAGGAATGA * NM:i:1 MD:Z:38 AS:i:31 XS:i:29 SA:Z:gi|224384768|gb|CM000663.1|,3668552,-,24S29M1D9M,1,2; #INDEL_lower_path_17964|P_1:30_10_8|low|nb_pol_1|left_unitig_length_346|right_unitig_length_815|C1_20|C2_43|G1_0/1:321,17,162|G2_1/1:848,120,10|rank_0.46189 2064 gi|224384768|gb|CM000663.1| 3668552 1 24H29M1D9M * 0 0 TTTTTTTCTTTTTTTTCTTTTTTTATTTCTTTTTTCTT * NM:i:2 MD:Z:12C16^T9 AS:i:30 XS:i:26 SA:Z:gi|224384768|gb|CM000663.1|,191102966,+,23S8M1I30M,1,1; XA:Z:gi|224384768|gb|CM000663.1|,-197957308,27S26M9S,0; #--------------------------------------------------------------------------------------------------------------------------- def CheckCoupleVariantID(self): """Test if the couple of paths has the same ID""" IDVariantUp = self.upper_path.listSam[0].split("_")[3] IDVariantLow= self.lower_path.listSam[0].split("_")[3] bitwiseFlag=int(self.upper_path.listSam[1]) if IDVariantUp != IDVariantLow: if bitwiseFlag >0: #& 2048 : #Checks if it's a supplementary alignment print("Supplementary alignment:") print(self.upper_path.listSam) return (2) else : print("WARNING two consecutive lines do not store the same variant id: ") print(self.upper_path.listSam) print(self.upper_path.listSam[1]) print(self.lower_path.listSam) return (2) else: return (0) ############################################################################################# ############################################################################################# class PATH(): """corresponds to one path of a discoSnp prediction""" def __init__(self,line): self.listCoverage=[] #list of all the coverage by sample for the path self.dicoMappingPos={} #dictionnary with all the mapping positions associated with their number of mismatches with the reference self.listNucleotideReverse=[] #list of all the variant (snp) of the path on the reverse strand self.listNucleotideForward=[] #list of all the variant (snp) of the path on the forward strand self.boolReverse=None #Boolean to know if the strand is reverse(-1) or forward(1) self.posMut=None #MD tag of the samfile "MD:Z:5A10A0A25G17" => 5A10A0A25G17 self.cigarcode=None #cigarcode of the samfile "61M" self.boolRef=None #Boolean to know if the path is identical to the reference self.nucleoRef=None #Nucleotide corresponding to the variant on the reference self.nucleo=None #nucleotide corresponding of the variant on the path self.listPosVariantOnPathToKeep=[] #list of the positions of all the variant on the in case of Reverse or Forward mapped path self.listPosReverse=[] #mapping position(s) of the variant on the path (if it is mapped on the reverse strand) self.listPosForward=[] #mapping position(s) of the variant on the path (if it is mapped on the forward strand) self.correctedPos=0 #list or position of the mapping variant by taking into account the shift with the reference self.listFQQuality=[] #string of all the quality scores of every variant if ">" not in line: # Case of samfile self.listSam=line.rstrip('\r').rstrip('\n').split('\t') self.discoName = self.listSam[0] self.seq = self.listSam[9] #gets the sequence of the path self.mappingPosition=abs(int(self.listSam[3])) #mapping position of the path self.RetrieveDicoMappingPosition() self.CheckBitwiseFlag() else: #Mode ghost fastafile line=line.strip('>').split("\n") line.pop() self.listSam=line #gets the sequence of the path self.discoName = self.listSam[0] self.seq="" self.mappingPosition=0 #mapping position of the path self.boolReverse="." #We need to define the absence of strand def RetrieveXA(self,VCFObject): localXA="" for position,(_,cigarcode) in self.dicoMappingPos.items(): if cigarcode!="": chromosome = '_'.join(position.split('_')[:-1]) mapping_position = int(position.split('_')[-1]) shifted_mapping_position = shift_from_cigar_code(cigarcode,abs(mapping_position)+self.listPosVariantOnPathToKeep[0]-1) localXA+=chromosome+'_'+str(shifted_mapping_position)+"," localXA=localXA.rstrip(',') VCFObject.XA = localXA #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def RetrieveSeq(self,seq): """Getter for sequence: fills path object""" self.seq=seq def RetrieveDicoMappingPosition(self): """Retrieves for each path alignment information in a list ; retrieves a dictionary with all the positions of a path and the number of associated mismatch""" variant = self.listSam listXA = None strXA = None alternative_positions = None nbMismatch = None variant_position = abs(int(variant[3])) variant_chromosome = variant[2] #Error list with mapping positions very close to the first position given by bwa # listerreur=set([(int(variant[3])-1),(int(variant[3])+1),(int(variant[3])+2),(int(variant[3])+3),(int(variant[3])-3),(int(variant[3])-2),int(variant[3])]) #Creation of a dict with mapping position associated with number of mismatch # XA cannot occur in the 11 mandatory fields: QNAME FLAG RNAM POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL if 'XA:Z' in ''.join(variant[11:]): # XA: tag for multiple mapping : Checks if the upper path is multiple mapped : XA Alternative hits; format: (chr,pos,CIGAR,NM;)* for item in variant[11:]: if "XA:Z" in item: #Parsing XA tag listXA=item.split(":")[2].split(';') strXA = ','.join(listXA) # TODO: useless (?) listXA = strXA.split(',') # TODO: useless (?) listXA.pop() alternative_positions = listXA #position=[chrom1,pos1,cigarcode1,number of mismatch1 , chrom2,pos2,cigarcode2,number of mismatch2,...]. pos_i may be negative, in case of revcomp mapping. self.XA = item break #no need to search for XA in other fields i=1 while i < len(alternative_positions): #Runs through the list 4 by 4 to get all the positions if alternative_positions[i-1]==variant_chromosome and\ abs(int(alternative_positions[i])) > variant_position-4 and\ abs(int(alternative_positions[i])) < variant_position+4: i+=4 continue #Checks if the position is not too close to the main one self.dicoMappingPos[alternative_positions[i-1]+"_"+alternative_positions[i]]=[int(alternative_positions[i+2]), alternative_positions[i+1]]#the position is associated to the number of mismatch in a dictionary i+=4 if variant_position>0:#adds the main mapping position to the dictionary of all mapping positions posMut,nbMismatch = self.GetTag() #In case of mapped variant without MD TAG : self.dicoMappingPos[abs(int(variant[3]))]=[nbMismatch,""] #--------------------------------------------------------------------------------------------------------------------------- # #FLAG = field to test #1 read paired #2 read mapped in proper pair #4 read unmapped #8 mate unmapped #16 read reverse strand #32 mate reverse strand #64 first in pair #128 second in pair #256 not primary alignment #512 read fails platform/vendor quality checks #1024 read is PCR or optical duplicate #2048 supplementary alignment # #--------------------------------------------------------------------------------------------------------------------------- def CheckBitwiseFlag(self): """Checks if the BitwiseFlag contains the tested value such as : read reverse strand, read unmmaped and so on.""" if int(self.listSam[1]) & 16:#Reverse strand self.boolReverse="-1" self.listPosVariantOnPathToKeep = self.listPosReverse return if int(self.listSam[1]) & 4: #Unmapped self.listPosVariantOnPathToKeep = self.listPosForward self.boolReverse="." return #else: #Forward strand self.listPosVariantOnPathToKeep = self.listPosForward self.boolReverse="1" #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def CigarcodeChecker(self): """Checks in the cigarcode of the samfile if there is a shift in the alignment between the path and the reference""" cigarcode = self.listSam[5] parsingCigarCode=re.findall(r'(\d+|[A-Za-z])',cigarcode) #ParsingCigarCode=['2', 'S', '3', 'M', '1', 'I', '25', 'M'] listPosRef=[] listShift=[] somme=0 shift=0 pos=0 i=1 j=0 listpol = self.listPosVariantOnPathToKeep while i=int(listpol[0]): posRef=int(listpol[0])+shift#takes into account the shift to add the after the mapping position and the variant position in the sequence listShift.append(shift) listPosRef.append(posRef) return(listPosRef,listShift) elif len(listpol)>1:#Close SNPs while int(pos)>=int(listpol[j]):#Goes through the list of position (close snps) and see if the position is affected by the shift (means shift before the position)self.listNucleotideForward posRef=int(listpol[j])+shift #Add the shift to the position (to get the real position of the snp on the reference) listPosRef.append(posRef) #Add the position to the list by taking into account the shift only if the current position listShift.append(shift) if j<(len(listpol)-1): j+=1 else: return(listPosRef,listShift) i+=2 #--------------------------------------------------------------------------------------------------------------------------- #Example of unmmaped variant : soft clipping of the variant #['SNP_higher_path_146392|P_1:30_A/G,P_2:45_C/T,P_3:48_T/G|high|nb_pol_3|left_unitig_length_157|right_unitig_length_564|C1_13|C2_1|G1_0/1:389,20,170|G2_1/1:828,117,10|rank_0.43053', '0', 'gi|224384768|gb|CM000663.1|', '229146041', '52', '79M', '*', '0', '0', 'CTTTCTATCTCAAAAGCAGCCACAGACCACATGTAAACAAATAAGCGGTGCCATGTTCCAATAAAACTTTATTTACAGA', '*', 'NM:i:5', 'MD:Z:15T23G1G4T23G8', 'AS:i:54', 'XS:i:30'] #['SNP_lower_path_146392|P_1:30_A/G,P_2:45_C/T,P_3:48_T/G|high|nb_pol_3|left_unitig_length_157|right_unitig_length_564|C1_24|C2_42|G1_0/1:389,20,170|G2_1/1:828,117,10|rank_0.43053', '16', 'gi|224384768|gb|CM000663.1|', '15779841', '25', '30M49S', '*', '0', '0', 'TCTGTAAATAAAGTTTTATTGGAACATGGCCCCACTTATTTGTTTACACGTGGTCTGTGGCTGCTTTTGAGATAGAAAG', '*', 'NM:i:0', 'MD:Z:30', 'AS:i:30', 'XS:i:25', 'XA:Z:gi|224384768|gb|CM000663.1|,+7370519,52S27M,1;gi|224384768|gb|CM000663.1|,-95561814,27M52S,1;'] #SNP_higher_path_215581|P_1:30_C/G|low|nb_pol_1|left_unitig_length_8|right_unitig_length_1|C1_7|C2_17|G1_0/1:554,48,75|G2_0/1:268,13,248|rank_0.32064 0 gi|224384768|gb|CM000663.1| 232979913 7 31S30M * 0 0 TCAAGACCAGCCTAGGCAACATAGAGATACCATGTCTCTACAAAAAATTAAAAAAAAAAAA * NM:i:0 MD:Z:30 AS:i:30 XS:i:28 XA:Z:gi|224384768|gb|CM000663.1|,-241321283,29M32S,1;gi|224384768|gb|CM000663.1|,-114672467,28M33S,0; #SNP_lower_path_215581|P_1:30_C/G|low|nb_pol_1|left_unitig_length_8|right_unitig_length_1|C1_31|C2_18|G1_0/1:554,48,75|G2_0/1:268,13,248|rank_0.32064 16 gi|224384768|gb|CM000663.1| 65214684 37 59M2S * 0 0 TTTTTTTTTTTTAATTTTTTGTAGAGACATCGTATCTCTATGTTGCCTAGGCTGGTCTTGA * NM:i:3 MD:Z:16A23A5A12 AS:i:44 XS:i:30 #--------------------------------------------------------------------------------------------------------------------------- def ReferenceChecker(self,shift,posCentraleRef,VCFObject,PosVariant): """Function which allows to get the MD tag parsing; checks if path nucleotide is identical to the reference nucleotide""" i=0 posMut,nbMismatch = self.GetTag() boolDel=True pos=shift #Allows or disallows soft clip in BWA #if int(shift)<=-(int(PosVariant)):#Test if the variant is really mapped (soft clipping > variant position => unmmaped variant) if int(shift)<=-2:#we allow 2 soft clip (for bwa mem) self.boolRef=False self.mappingPosition=0#It is considered that the path is unmapped return() nucleoRef=None matchInt=None parsingPosMut=re.findall(r'(\d+|[A-Za-z]|\^)',posMut) while i it means that the nucleotide is different in the variant and in the reference self.boolRef=False self.nucleoRef=parsingPosMut[i] break else: #If the last item of the list of the MD tag is an integer => it means that the nucleotide of the allele is identical to the reference self.boolRef=True break if pos>posCentraleRef: #If the current position is bigger than the variant position it means that the nucleotide of the variant is identical to the reference self.boolRef=True break i+=1 if pos0:#Check if the variant is really mapped if "MD:" not in str(variant):#Not MD Tag in the variant we deduce the value from the cigarcode print ("!!! No MD tag in your sam file : Could you try with the last version of bwa (upper than 0.7.8) ?") sys.exit() else: for field in self.listSam[11:]: if "MD:" in field: posMut = field.split(":")[2] #MD tag parsing break return (posMut,nbMismatch) #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def CheckPosVariantFromRef(self,VCFObject): #Validation snp first part """Checks if the variant is identical to the reference or not ; defines the nucleotide on the reference""" dicoClose={} nucleo=None boolRef=None nucleoRef=None if int(self.mappingPosition)>0: #Gets the shift by positions (insertion,deletion,sofclipping) and update of the position on the path listCorrectedPos,listShift = self.CigarcodeChecker() self.correctedPos=listCorrectedPos #Defines if the path is identical to the reference and what is the nucleotide on the reference i=0 for i in range(len(listCorrectedPos)):#Loops on the list of corrected positions self.ReferenceChecker(listShift[i],listCorrectedPos[i],VCFObject,self.listPosVariantOnPathToKeep[i])#Checks if the path is identical to the reference genome if int(self.mappingPosition)<=0:# Case => variant considered as unmapped because of soft clipping so we have to check again if the mapping position break if self.boolReverse=="1" and self.listNucleotideForward!=[]:#If we are on the forward strand => defines the nucleotide for the current snp or indel. self.nucleo = self.listNucleotideForward[i] if self.nucleoRef==None:#If there is no reference nucleotide given by ReferenceChecker, it means that the variant is equal to the reference so we defined it ! self.nucleoRef = self.listNucleotideForward[i] elif self.boolReverse=="-1" and self.listNucleotideReverse!=[]:#If we are on the reverse strand => defines the nucleotide for the current snp or indel. self.nucleo = self.listNucleotideReverse[i] if self.nucleoRef==None:#If there is no reference nucleotide given by ReferenceChecker, it means that the variant is equal to the reference so we defined it ! self.nucleoRef = self.listNucleotideReverse[i] dicoClose[self.listPosVariantOnPathToKeep[i]]=[(self.boolRef),(self.nucleoRef),(listCorrectedPos[i]),(self.nucleo),(self.boolReverse),(int(listCorrectedPos[i])+int(self.mappingPosition))]#Fills the dictionary to keep all informations for close snps if i==0: #Keeps the information of the first snp/indel to fill the attributs of the path object boolRef = self.boolRef nucleoRef = self.nucleoRef nucleo = self.nucleo self.nucleoRef=None#Resets the reference nucleotide for the next snp (case of close snps) if int(self.mappingPosition)<=0:#Case of unmapped path listCorrectedPos = self.listPosForward#We keep the forward position for every snps/indel on the path self.listPosVariantOnPathToKeep = self.listPosForward shift=0#There is no shift with the reference self.boolReverse="." self.nucleoRef="." self.boolRef=False boolRef=False nucleoRef='.' self.correctedPos=listCorrectedPos if self.listNucleotideForward!=[]: self.nucleo = self.listNucleotideForward[0] i=0 for i in range(len(self.listNucleotideForward)):#Loops on the list of position to fill the dictionary for close snps dicoClose[self.listPosForward[i]]=[(False),(self.nucleoRef),(listCorrectedPos[i]),(self.listNucleotideForward[i]),self.boolReverse,(int(listCorrectedPos[i])+int(self.mappingPosition))] if i==0: boolRef=False nucleoRef = self.nucleoRef nucleo = self.nucleo #Fills the attribut of the path with the information of the first snp/indel self.boolRef=boolRef self.nucleoRef=nucleoRef self.nucleo=nucleo return(dicoClose) ############################################################################################# ############################################################################################# class SNP(VARIANT): def __init__(self,line1,line2): VARIANT.__init__(self,line1,line2) #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def WhichPathIsTheRef(self,VCFObject): """Finds which path is identical to the reference genome (with boolRef) and defines it as the ref : specific method for each type of variant""" VARIANT.WhichPathIsTheRef(self,VCFObject) posUnmapped = max(self.len_unitig_left,self.len_contig_left) #Takes into account the length of the unitig/contig for the position of unmapped allele (position of the allele on the lower path) # self.CheckContigUnitig(self.len_unitig_left,self.len_contig_left) #Takes into account the length of the unitig/contig for the position of unmapped allele (position of the allele on the lower path) #--------------------------------------------------------------------------------------------------------------------------- ##Case : two mapped paths if self.upper_path.mappingPosition>0 and self.lower_path.mappingPosition>0: ##The path identical to the reference is the lower path if self.lower_path.boolRef == True and self.upper_path.boolRef==False: VCFObject.chrom = self.lower_path.listSam[2] VCFObject.ref = self.lower_path.nucleo VCFObject.alt = self.upper_path.nucleo VCFObject.reverse = self.lower_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef ##The path identical to the reference is the upper path elif self.upper_path.boolRef == True and self.lower_path.boolRef==False: VCFObject.chrom = self.upper_path.listSam[2] VCFObject.ref = self.upper_path.nucleo VCFObject.alt = self.lower_path.nucleo VCFObject.reverse = self.upper_path.boolReverse VCFObject.nucleoRef = self.upper_path.nucleoRef ##No path is identical to the reference => lexicographique choice elif self.upper_path.boolRef == False and self.lower_path.boolRef==False: if self.upper_path.nucleo < self.lower_path.nucleo: VCFObject.chrom = self.upper_path.listSam[2] VCFObject.ref = self.upper_path.nucleo VCFObject.alt = self.lower_path.nucleo VCFObject.reverse = self.upper_path.boolReverse VCFObject.nucleoRef = self.upper_path.nucleoRef elif self.upper_path.nucleo > self.lower_path.nucleo: VCFObject.chrom = self.lower_path.listSam[2] VCFObject.ref = self.lower_path.nucleo VCFObject.alt = self.upper_path.nucleo VCFObject.reverse = self.lower_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef #--------------------------------------------------------------------------------------------------------------------------- ##Case : Lower path mapped and upper path unmapped elif self.upper_path.mappingPosition<=0 and self.lower_path.mappingPosition>0: VCFObject.chrom = self.lower_path.listSam[2] VCFObject.ref = self.lower_path.nucleo VCFObject.alt = self.upper_path.nucleo VCFObject.reverse = self.lower_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef #--------------------------------------------------------------------------------------------------------------------------- ##Case : Upper path mapped and lower path unmapped elif self.upper_path.mappingPosition>0 and self.lower_path.mappingPosition<=0: VCFObject.chrom = self.upper_path.listSam[2] VCFObject.ref = self.upper_path.nucleo VCFObject.alt = self.lower_path.nucleo VCFObject.reverse = self.upper_path.boolReverse VCFObject.nucleoRef = self.upper_path.nucleoRef #--------------------------------------------------------------------------------------------------------------------------- ##Case : Both paths are unmapped elif self.upper_path.mappingPosition<=0 and self.lower_path.mappingPosition<=0: self.mappingPositionCouple=int(posUnmapped) if self.lower_path.nucleolen(self.upper_path.seq): #if the sequence of the upper path is the smallest self.lower_path.nucleo = self.insertReverse self.upper_path.nucleo = self.ntStartReverse else: self.upper_path.nucleo = self.insertReverse self.lower_path.nucleo = self.ntStartReverse #If the upper path or the lower path is the ref and if it is a reference mappind # 7/6/2016: CL and PP: changed this elif for the else. Avoids an issue with self.upper_path.nucleo not initialized. Large tests results are the same. else:#if (self.lower_path.boolRef==True and (self.lower_path.boolReverse=="1" or self.lower_path.boolReverse==".")) or (self.upper_path.boolRef==True and (self.upper_path.boolReverse=="1" or self.lower_path.boolReverse==".")): if len(self.lower_path.seq)>len(self.upper_path.seq): #if the sequence of the upper path is the smallest : self.lower_path.nucleo = self.insertForward self.upper_path.nucleo = self.ntStartForward else: self.upper_path.nucleo = self.insertForward self.lower_path.nucleo = self.ntStartForward #In case of mapped variant if old_boolRef_Low==None and old_boolRef_Up==None: #In case of unmapped variant or too much soft clip we have an other way to fill vcf ##Fills the VCF if the upper path is considered as the reference if self.upper_path.boolRef==True: if len(self.upper_path.nucleo) == len(self.insertForward): self.upper_path.nucleoRef="." VCFObject.variantType="DEL" else: self.upper_path.nucleoRef="." VCFObject.variantType="INS" if self.upper_path.mappingPosition>0: VCFObject.chrom = self.upper_path.listSam[2] else: VCFObject.chrom = self.upper_path.listSam[0].split("|")[0] VCFObject.ref = self.upper_path.nucleo VCFObject.alt = self.lower_path.nucleo VCFObject.reverse = self.upper_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef ##Fills the VCF if the lower path is considered as the reference elif self.lower_path.boolRef==True: if len(self.lower_path.nucleo)==len(self.insertForward): self.lower_path.nucleoRef="." VCFObject.variantType="DEL" else: self.lower_path.nucleoRef="." VCFObject.variantType="INS" if self.lower_path.mappingPosition>0: VCFObject.chrom = self.lower_path.listSam[2] else: VCFObject.chrom = self.lower_path.listSam[0].split("|")[0] VCFObject.ref = self.lower_path.nucleo VCFObject.alt = self.upper_path.nucleo VCFObject.reverse = self.lower_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef # Unmapped variants or too much soft clipped (=> the variant will be considered as unmapped) if (self.upper_path.mappingPosition<=0 and self.lower_path.mappingPosition<=0) or (old_boolRef_Low==False and old_boolRef_Up==False): if self.lower_path.boolRef==True: self.mappingPositionCouple=int(posUnmapped) else: self.mappingPositionCouple = int(posUnmapped) if len(self.upper_path.nucleo) == len(self.insertForward): VCFObject.chrom = self.upper_path.discoName.split("|")[0] self.upper_path.boolRef = True self.upper_path.nucleoRef = "." VCFObject.variantType = "DEL" VCFObject.ref = self.upper_path.nucleo VCFObject.alt = self.lower_path.nucleo VCFObject.reverse = self.upper_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef else: VCFObject.chrom = self.lower_path.discoName.split("|")[0] self.upper_path.boolRef = False self.upper_path.nucleoRef = "." VCFObject.variantType = "INS" VCFObject.ref = self.lower_path.nucleo VCFObject.alt = self.upper_path.nucleo VCFObject.reverse = self.lower_path.boolReverse VCFObject.nucleoRef = self.lower_path.nucleoRef def FillVCF(self,VCFfile,nbGeno,table,VCFObject): table[4]=str(VCFObject.alt) VARIANT.FillVCF(self,VCFfile,nbGeno,table,VCFObject) #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def RetrieveMappingPositionCouple(self): VARIANT.RetrieveMappingPositionCouple(self) self.mappingPositionCouple=int(self.mappingPositionCouple) ############################################################################################# ############################################################################################# class SNPSCLOSE(VARIANT): def __init__(self,line1,line2): VARIANT.__init__(self,line1,line2) self.dicoCloseSNPUp={} #dictionnary with all the informations for close snps : boolean to know if the allele is identical to the reference,position on the variant, reverse nucleotide for the allele, if the path is reverse or not, mapping position self.dicoCloseSNPLow={} #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def RetrieveDicoClose(self,dicoCloseUp,dicoCloseLow):# Gets the outputs of the method CheckPosVariantFromRef (PATH CLASS) self.dicoCloseSNPUp=dicoCloseUp self.dicoCloseSNPLow=dicoCloseLow #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def WhichPathIsTheRef(self,VCFObject): VARIANT.WhichPathIsTheRef(self,VCFObject) VCFObject.phased=True table = [0] * 10 #Create a 10 cols array tablebis=[] posUnmapped = max(self.len_unitig_left,self.len_contig_left) listPositionPolymorphismeOnPathUp = self.upper_path.listPosVariantOnPathToKeep listPositionPolymorphismeOnPathLow = self.lower_path.listPosVariantOnPathToKeep VCFObject.nucleoRef=[] listPolymorphismPos=[] ##Case : two mapped paths if self.upper_path.mappingPosition>0 and self.lower_path.mappingPosition>0: ##Sorts the list of position to get the smallest one and its position in the list unsorted! ##Keeps the close snps to sort them : indeed all the lists and dictionnaries :listnucleoUp,listPositionPolymorphismeOnPathUp,listPositionPolymorphismeOnPathLow,listnucleoLow self.dicoCloseSNPUp,self.dicoCloseSNPLow are classified according to dicoHeader so if we start by sorting we lose the correspondence between data listSortedPosUp=list(listPositionPolymorphismeOnPathUp) listSortedPosUp.sort() indexSmallestPosUp=listPositionPolymorphismeOnPathUp.index(listSortedPosUp[0]) listSortedPosLow=list(listPositionPolymorphismeOnPathLow) listSortedPosLow.sort() indexSmallestPosLow=listPositionPolymorphismeOnPathLow.index(listSortedPosLow[0]) self.upper_path.boolRef = self.dicoCloseSNPUp[listPositionPolymorphismeOnPathUp[indexSmallestPosUp]][0] self.lower_path.boolRef = self.dicoCloseSNPLow[listPositionPolymorphismeOnPathLow[indexSmallestPosLow]][0] #self.upper_path.boolRef = self.dicoCloseSNPUp[listPositionPolymorphismeOnPathUp[indexSmallestPosUp]][0] #self.lower_path.boolRef = self.dicoCloseSNPLow[listPositionPolymorphismeOnPathLow[indexSmallestPosLow]][0] #Decides what is the smallest position according to the reference path (useful if the paths are not aligned on the same strand) if (self.upper_path.boolRef==False and self.lower_path.boolRef==False) or (self.upper_path.boolRef==True and self.lower_path.boolRef==True): VARIANT.MismatchChecker(self) if self.upper_path.boolRef==True and self.lower_path.boolRef==False: #The smallest position identical to the reference is on the upper path listPolymorphismPos=listPositionPolymorphismeOnPathUp elif self.lower_path.boolRef==True and self.upper_path.boolRef==False:#The smallest position identical to the reference is on the lower path listPolymorphismPos=listPositionPolymorphismeOnPathLow elif self.upper_path.boolRef==False and self.lower_path.boolRef==False: #Both paths are different from the reference => choice with the lexicographical order if self.upper_path.nucleo0 and self.lower_path.mappingPosition<=0: comptPol=0 VCFObject.chrom = self.upper_path.listSam[2] VCFObject.reverse = self.upper_path.boolReverse for comptPol in range(len(listPositionPolymorphismeOnPathUp)): if (int(self.upper_path.boolReverse)==-1): nucleoLow=ReverseComplement(self.lower_path.listNucleotideForward[comptPol]) nucleoUp = self.upper_path.listNucleotideReverse[comptPol] elif int(self.upper_path.boolReverse)==1: nucleoLow = self.lower_path.listNucleotideForward[comptPol] nucleoUp = self.upper_path.listNucleotideForward[comptPol] nucleoRefLow="." positionSnpUp=int(self.dicoCloseSNPUp[listPositionPolymorphismeOnPathUp[comptPol]][5])-1 nucleoRefUp = self.dicoCloseSNPUp[listPositionPolymorphismeOnPathUp[comptPol]][1] table[1]=positionSnpUp table[3]=nucleoUp table[4]=nucleoLow VCFObject.nucleoRef.append([nucleoRefUp,positionSnpUp]) tablebis.append(list(table)) tablebis=sorted(tablebis, key=lambda colonnes: colonnes[1]) return (tablebis) #--------------------------------------------------------------------------------------------------------------------------- ##Case : Lower path mapped and upper path unmapped elif self.upper_path.mappingPosition<=0 and self.lower_path.mappingPosition>0: VCFObject.chrom = self.lower_path.listSam[2] VCFObject.reverse = self.lower_path.boolReverse for comptPol in range(len(listPositionPolymorphismeOnPathLow)): if (int(self.lower_path.boolReverse)==-1): nucleoUp=ReverseComplement(self.upper_path.listNucleotideForward[comptPol]) nucleoLow = self.lower_path.listNucleotideReverse[comptPol] elif int(self.lower_path.boolReverse)==1: nucleoUp = self.upper_path.listNucleotideForward[comptPol] nucleoLow = self.lower_path.listNucleotideForward[comptPol] nucleoRefUp="." positionSnpLow=int(self.dicoCloseSNPLow[listPositionPolymorphismeOnPathLow[comptPol]][5])-1 nucleoRefLow = self.dicoCloseSNPLow[listPositionPolymorphismeOnPathLow[comptPol]][1] table[1] = positionSnpLow table[3] = nucleoLow table[4] = nucleoUp VCFObject.nucleoRef.append([nucleoRefLow,positionSnpLow]) tablebis.append(list(table)) tablebis=sorted(tablebis, key=lambda colonnes: colonnes[1]) return (tablebis) #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def FillVCF(self,VCFfile,nbGeno,table,VCFObject): """print the VCFFile for each line of close snps""" #ID=1 subIDs = range(1,len(table)+1,1) #subIDs for the different SNPs of a given close bubble : ex: 454_1, 454_2 and 454_3 for 3 SNPs inside the bubble of ID 454 if VCFObject.reverse == "-1": # if the bubble of close SNPs is mapped on reverse strand, their subIDs are reversed, this way the subID always correspond to the order of apparition of the SNP in the bubble path subIDs = range(len(table),0,-1) i=0 nucleoRef=None if VCFObject.chrom=="*": table[line][0] = self.listSam[0].split("_")[0] for line in range(len(table)): for i in range(len(VCFObject.nucleoRef)): if table[line][1] in VCFObject.nucleoRef[i]: nucleoRef=VCFObject.nucleoRef[i][0] table[line][0]=VCFObject.chrom table[line][2]=str(self.variantID)+"_"+str(subIDs[line]) table[line][5]="." table[line][6]=VCFObject.filterField table[line][7]=f"Ty={VCFObject.variantType};Rk={self.rank};UL={self.len_unitig_left};UR={self.len_unitig_right};CL={self.len_contig_left};CR={self.len_contig_right};Genome={nucleoRef};Sd={VCFObject.reverse}" if VCFObject.filterField=="MULTIPLE" and VCFObject.XA: table[line][7]+=f";XA={VCFObject.XA}" #TODO: eviter ces "replace" # table[line][7]=table[line][7].replace("None",".") # done during PrintOneLine table[line][7]=table[line][7].replace("none",".") table[line][7]=table[line][7].replace("=;","=.;") table[line][8]=VCFObject.formatField table[line][9]=VCFObject.genotypes #ID+=1 i+=1 error=VCFObject.CheckOutputConsistency(table,self) if error == 0: for l in range(len(table)): VCFObject.PrintOneLine(table[l],VCFfile) ############################################################################################# ############################################################################################# class VCFFIELD(): def __init__(self): ##VCF Fields self.filterField=None self.variantType=None self.phased=None #phased genotype self.formatField=None self.genotypes="" #genotype field =>put it in the vcf self.chrom="" self.ref="" self.alt="" self.qual="" self.nucleoRef=None self.reverse="" self.XA=None #--------------------------------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------------------------------- def PrintOneLine(self,table,VCF): """Prints the line of the current SNP in the VCF file.""" if table==[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]: return() for i in range(len(table)): element=table[i] element=str(element).replace("None",".") # TODO: peut on éviter ce nouveau 'replace'... ? VCF.write((str(element)).strip()) if i0: sys.stderr.write(" !!! Line where the error occurred !!!\n") sys.stderr.write(str(VARIANT.upper_path.listSam)+"\n") sys.stderr.write(str(VARIANT.lower_path.listSam)+"\n") else : return (error)