#!/usr/bin/python ######################################################################################### ## Author: Pravech Ajawatanawong ## ## Last Update: 8 JUNE 2021 ## ## Program: SeqFIRE: Sequence Feature and Indel Region Extractor ## ## Version: 2.21 (2021) ## ######################################################################################### import getopt, sys, os, re, warnings from math import log10 #################################### ## D E F A U L T V A L U E S ## #################################### amino_acid = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T',\ 'V', 'W', 'X', 'Y', '-'] pMatrices = {'PAM250': [('?'), ('M','I','V','L'), ('D','N','H','Q','E'), ('F','I','L'), ('S','P','A'),\ ('S','A','G'), ('Q','K','N'), ('R','H','Q'), ('S','T'), ('R','K','Q'), ('S','N'),\ ('F','W'), ('R','W'), ('G','D')], 'PAM60': [('?'), ('S','A','T'), ('H','N'), ('S','N'),\ ('N','D'), ('D','E'), ('H','Q'), ('Q','E'), ('Y','F'), ('R','K'), ('M','I'),\ ('L','M'), ('I','V')], 'BLOSUM40': [('?'), ('S','T'), ('S','A'), ('S','Q','N'),\ ('H','Y','M'), ('D','E'), ('N','D'), ('H','N'), ('W','Y','F'), ('E','Q','K'),\ ('K','R'), ('R','Q'), ('L','M','I','V'), ('E','K'), ('A','G'), ('L','F','I'),\ ('V','T')], 'BLOSUM62': [('?'), ('S','T'), ('S','A'), ('S','N'), ('H','Y'), ('D','E'),\ ('N','D'), ('H','N'), ('W','Y','F'), ('E','Q','K'), ('K','R','Q'),\ ('L','M','I','V')], 'BLOSUM80': [('?'), ('Q','R','K'), ('Q','E','K'), ('E','D'), ('D','N'),\ ('Q','H'), ('Y','H'), ('Y','W'), ('Y','F'), ('S','T'), ('S','A'), ('M','V','L','I')]} prog_title = '# SeqFIRE: Sequence Feature and Indel Region Extractor\n# version 1.0.1 (c) 2011' infile = 'infile' analysis_mode = 1 similarity_threshold = [[75.0,100.0]] percent_similarity = [[75.0,100.0]] percent_accept_gap = 40.0 p_matrix = 'NONE' p_matrix_2 = 'NONE' inter_indels = 3 partial = 'True' blocks = 3 strick_combination = 'False' combine_with_indel = 'False' fuse = 4 multidata = 1 output_mode = 2 seqType = 'Protein' def usage(): text=""" SeqFIRE: Sequence Feature and Indl Region Extractor Version 1.0 Copyright 2011 URL: www.seqfire.org --------------------------------------------------- [ HELP MESSAGE ] for quick help type >>> python seqfire.py -h for quick started just type >>> python seqfire.py -i {inputfile with path} {other options} +-------------------+ | Options available | +-------------------+ Option Description -h display SeqFIRE help -i path and name of your infile -a analysis module 1 = indel region module 2 = conserved block module Default = 1 -c similarity threshold for indel regions identification (0.0 - 100.0) Default = 75.0 -d similarity threshold for conserved blocks identification (0.0 - 100.0) Default = 75.0 -j percent accept gap (0.0 - 100.0) Default = 40.0 -g substitution group for indel identification (PAM250, PAM60, BLOSUM40, BLOSUM62, BLOSUM80, or NONE) Default = NONE -k substitution group for cnserved block identification (PAM250, PAM60, BLOSUM40, BLOSUM62, BLOSUM80, or NONE) Default = NONE -b interindel space (1-10) Default = 3 -t twilight treatment (True or False) Default = True -p partial sequence treatment (True or False) Default = False -s space between 2 conserved blocks (1, 2, 3, ...) Default = 3 -f minimum size for non-conserved block (1, 2, 3, ...) Default = 4 -r strick combination (True or False) Default = False -e get indel matrix together with conserved region (True or False) Default = False -m multiple data analysis 1 = single dataset 2 = multiple datasets (batch mode) Default = 2 -o output mode 1 = show output on screen 2 = save output in a file 3 = both show output on screen & save in a file For more details, please read the SeqFIRE manual (www.seqfire.org/help). """ print (text) def parseFasta(records): rec = [] seqs = records.split('>') del seqs[0] for seq in seqs: if '\r\n' in records: p = seq.split('\r\n') else: p = seq.split('\n') if len(p[0]) > 60: title = p[0][:60] else: title = p[0] del p[0] seq = ''.join("%s" % (k) for k in p) rec.append([title, seq.upper()]) title_length = 0 for r in rec: if len(r[0]) > title_length: title_length = len(r[0]) if title_length < 23: for r in rec: r[0] = str(r[0]) + ' '*(23 - len(r[0])) + ': ' else: for r in rec: r[0] = str(r[0]) + ' '*(title_length - len(r[0])) + ': ' return rec def getGapProfile(seq_lists): gap_profile = '' for i in range(len(seq_lists[0][1])): g = [] for seq_list in seq_lists: g.append(seq_list[1][i]) if '-' in g: gap_profile = gap_profile + '-' else: gap_profile = gap_profile + 'X' #print("Gap profiling complete:\n",gap_profile) return gap_profile ######################################################################### ###### ###### ###### I N D E L R E G I O N M O D U L E ###### ###### ###### ######################################################################### def getBodyAlignmentIndexes(seq_lists): gap_profile = getGapProfile(seq_lists) ind_n = gap_profile.find('XXXXX') # Search for N-terminal edge of alignment body if ind_n==-1: ind_n=0 ind_reverse_c = gap_profile[::-1].find('XXXXX') # Search for C-terminal edge of alignment body ind_c = len(gap_profile) - ind_reverse_c - 1 return [ind_n, ind_c] def getBodyOfAlignment(ind, seq_lists): body_alignment = [] for seq_list in seq_lists: record = [] title = seq_list[0] seq = seq_list[1][ind[0]:ind[1]+1] record = [title, seq] body_alignment.append(record) return body_alignment def genMaskRef(seq_lists, position): ref = '' for i in range(0, position, 1): count = 0 for seq_list in seq_lists: if not seq_list[1][i] == '-': count += 1 if float(count)/float(len(seq_lists))*100.0 >= 55.0: ref = ref + '?' else: ref = ref + '-' return ref def getNRaggedRegion(n_position, seq_lists): n_ragged = [] if partial == 'True': if n_position <= 3: for seq_list in seq_lists: seq = '?' * n_position n_ragged.append([seq_list[0], seq]) else: maskRef = genMaskRef(seq_lists, n_position) for seq_list in seq_lists: rec = []; newSeq = [] for i in range(len(maskRef)): if maskRef[i]=='?': if seq_list[1][i]=='-': newSeq.append('?') else: newSeq.append(seq_list[1][i]) else: newSeq.append(seq_list[1][i]) rec = [seq_list[0], "".join(newSeq)] n_ragged.append(rec) else: for seq_list in seq_lists: seq = seq_list[1][:n_position] n_ragged.append([seq_list[0], seq]) return n_ragged def getCRaggedRegion(c_position, seq_lists): c_ragged = [] if partial == 'True': if len(seq_lists[0][1])-c_position-1 <= 3: for seq_list in seq_lists: seq = '?' * (len(seq_list[1]) - c_position - 1) c_ragged.append([seq_list[0], seq]) else: new_seq_lists = [] for seq_list in seq_lists: oneRecord = [seq_list[0], seq_list[1][::-1]] new_seq_lists.append(oneRecord) inv_maskRef = genMaskRef(new_seq_lists, len(new_seq_lists[0][1])-c_position-1) for seq_list in seq_lists: rec = []; newSeq = [] for i in range(len(new_seq_lists[0][1]) - c_position - 1): if inv_maskRef[i]=='?': if seq_list[1][i] == '-': newSeq.append('?') else: newSeq.append(seq_list[1][i]) else: newSeq.append(seq_list[1][i]) rec = [seq_list[0], ("".join(newSeq)[::-1])] c_ragged.append(rec) else: for seq_list in seq_lists: seq = seq_list[1][c_position+1:] c_ragged.append([seq_list[0], seq]) return c_ragged def genPseudoalignment(seq_lists): p, pseudoalignment = [], [] ragged_positions = getBodyAlignmentIndexes(seq_lists) body_of_alignment = getBodyOfAlignment(ragged_positions, seq_lists) # GENERATES THE N-TERMINAL RAGGED n_ragged_alignment = getNRaggedRegion(ragged_positions[0], seq_lists) # GENERATES THE C-TERMINAL RAGGED c_ragged_alignment = getCRaggedRegion(ragged_positions[1], seq_lists) # COMBINATION OF N-TERMINAL RAGGED, BODY OF ALIGNMENT AND C-TERMINAL RAGGED for k in range(0, len(body_of_alignment), 1): seq = [body_of_alignment[k][0], n_ragged_alignment[k][1] + body_of_alignment[k][1] + c_ragged_alignment[k][1]] pseudoalignment.append(seq) return pseudoalignment def getSimilarityScore(aa): if p_matrix == 'NONE': aa_dict = {'A': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0,\ 'K': 0, 'L': 0, 'M': 0, 'N': 0, 'P': 0, 'Q': 0, 'R': 0, 'S': 0,\ 'T': 0, 'V': 0, 'W': 0, 'X': 0, 'Y': 0, '?': 0} for a in aa: aa_dict[a] = aa_dict[a] + 1 k = max(aa_dict, key=aa_dict.get) similarity = (float(aa_dict[k]) * 100.0) / float(len(aa)) for onePair in similarity_threshold: if onePair[0]<=similarity<=onePair[1]: return 'C' return '.' else: pMatrix = pMatrices[p_matrix] c = [] aa_dict = {'A': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0,\ 'K': 0, 'L': 0, 'M': 0, 'N': 0, 'P': 0, 'Q': 0, 'R': 0, 'S': 0,\ 'T': 0, 'V': 0, 'W': 0, 'X': 0, 'Y': 0, '?': 0} for a in aa: aa_dict[a] = aa_dict[a] + 1 for aa_sets in pMatrix: m = 0 for aa_set in aa_sets: m = m + int(aa_dict[aa_set]) c.append((float(m)*100.0)/float(len(aa))) similarity = max(c) for onePair in similarity_threshold: if onePair[0]<=similarity<=onePair[1]: return 'C' return '.' def getConservedBlockProfile(pseudoSeq_lists): consBlkProfile = '' gapProfile = getGapProfile(pseudoSeq_lists) for i in range(0, len(gapProfile), 1): if gapProfile[i] == '-': consBlkProfile = consBlkProfile + '-' else: aa = [] for pseudoSeq_list in pseudoSeq_lists: aa.append(pseudoSeq_list[1][i]) ans = getSimilarityScore(aa) consBlkProfile = consBlkProfile + ans return consBlkProfile def getIndelProfile(conserved_block_profile): while (conserved_block_profile.count('-.')+conserved_block_profile.count('.-')>0): conserved_block_profile = conserved_block_profile.replace('-.', '--') conserved_block_profile = conserved_block_profile.replace('.-', '--') conserved_block_profile = conserved_block_profile.replace('.', 'C') for i in range(inter_indels-1, 0, -1): conserved_block_profile = conserved_block_profile.replace('-'+'C'*i+'-', '-'*(i+2)) m = re.search('C', conserved_block_profile) start_point = m.start() m = re.search('C', conserved_block_profile[::-1]) end_point = len(conserved_block_profile) - (m.end() - 1) sub_conserved_block_profile = conserved_block_profile[start_point:end_point] sub_conserved_block_profile = sub_conserved_block_profile.replace('-', 'I') sub_conserved_block_profile = sub_conserved_block_profile.replace('C', '.') return '.'*start_point + sub_conserved_block_profile + '.' * (len(conserved_block_profile) - end_point) def getIndelPositions(indel_profile): indel_position = [] for i in range(0, len(indel_profile)-1, 1): current_position = indel_profile[i] next_position = indel_profile[i+1] if current_position == '.' and next_position == 'I': indel_position.append(i+1) elif current_position == 'I' and next_position == '.': indel_position.append(i+1) return indel_position def getRuler(title_length, seq_length): scale_of_10 = '---------|' ruler = ' ' * title_length for i in range(50, seq_length, 50): ruler = ruler + ' ' * (50 - len(str(i))) + str(i) scale = ' ' * title_length + scale_of_10 * (seq_length//10) + scale_of_10[0:seq_length % 10] return ruler + '\n' + scale def genIndelAlignment(seq_lists, indel_profile, ruler, indel_positions): output_indel_1 = [] heading = '%s\n# OUTPUT: ANNOTATED ALIGNMENT\n# There are %s indels found.\n' % (prog_title, len(indel_positions)//2) output_indel_1.append(heading) output_indel_1.append(ruler) for seq_list in seq_lists: line = seq_list[0] + seq_list[1] output_indel_1.append(line) indelRef_line = ' '*len(seq_lists[0][0]) + indel_profile output_indel_1.append(indelRef_line) return output_indel_1 def genHomologAlignment(seq_lists, indel_profile, ruler, indel_positions): output_indel_4 = [] masked_alignment = [] homolog_positions = [] g = getBodyAlignmentIndexes(seq_lists) homolog_positions.append(g[0]) for j in indel_positions: homolog_positions.append(j) homolog_positions.append(g[1]) for seq_list in seq_lists: new_seq = '' for i in range(0, len(homolog_positions), 2): new_seq = new_seq + seq_list[1][homolog_positions[i]:homolog_positions[i+1]] masked_alignment.append([seq_list[0], new_seq]) heading = '''#NEXUS BEGIN DATA; \tDIMENSIONS NTAX=%s NCHAR=%s; \tFORMAT MISSING=? DATATYPE=PROTEIN GAP=-; \tOPTIONS GAPMODE=MISSING; MATRIX''' % (len(masked_alignment), len(masked_alignment[0][1])) output_indel_4.append(heading) #ruler = getRuler(len(masked_alignment[0][0]), len(masked_alignment[0][1])) r = getRuler(len(masked_alignment[0][0]), len(masked_alignment[0][1])).split('\n') ruler = '[' + r[0][1:] + ' '*(len(r[1])-len(r[0])) + ']\n[' + r[1][1:] + ']' output_indel_4.append(ruler) for s in masked_alignment: line = s[0] + s[1] output_indel_4.append(line) line = 'END;' output_indel_4.append(line) return output_indel_4 def search_for_simple_indels(pseudoalignments, inter_indels, indel_positions): indels = '' for i in range(0, len(indel_positions), 2): simple = True for pseudoalignment in pseudoalignments: indel_body = pseudoalignment[1][indel_positions[i]:indel_positions[i+1]] if not ('-'*len(indel_body) == indel_body or not '-' in indel_body): simple = False if simple: indels = indels + 's' else: indels = indels + 'c' return indels def genIndelLists(pseudoalignments, inter_indels, indel_positions, simple_indel_positions): output_indel_2 = [] heading = '%s\n# OUTPUT: INDEL LIST\n# There are %s indels found.\n# ** indicates masked sequence\n' % (prog_title, len(indel_positions)//2) output_indel_2.append(heading) id = 1 for i in range(0, len(indel_positions), 2): if simple_indel_positions[id-1] == 's': indel_type = 'simple indel' elif simple_indel_positions[id-1] == 'c': indel_type = 'complex indel' if indel_positions[i] + 1 == indel_positions[i+1]: line = '//\nIndel number: ' + str(id) + '\nIndel location in alignment: ' + str(indel_positions[i] + 1) + '\nSize of indel: ' +\ str(indel_positions[i+1] - indel_positions[i]) + ' alignment columns\nType: ' + indel_type + '\n' else: line = '//\nIndel number: ' + str(id) + '\nIndel location in alignment: ' + str(indel_positions[i] + 1) + '-' +\ str(indel_positions[i+1]) + '\nSize of indel: ' + str(indel_positions[i+1] - indel_positions[i]) +\ ' alignment columns\nType: ' + indel_type + '\n' for pseudoalignment in pseudoalignments: upstream_region = pseudoalignment[0] + pseudoalignment[1][indel_positions[i]-inter_indels:indel_positions[i]] indel_body = pseudoalignment[1][indel_positions[i]:indel_positions[i+1]] downstream_region = pseudoalignment[1][indel_positions[i+1]:indel_positions[i+1]+inter_indels] if '?' in upstream_region or '?' in indel_body or '?' in downstream_region: line = line + '\n' + upstream_region + ' ' + indel_body + ' ' + downstream_region + ' **' else: line = line + '\n' + upstream_region + ' ' + indel_body + ' ' + downstream_region output_indel_2.append(line) id += 1 return output_indel_2 def getIndelMatrix(output_indel_2, simple_indel_positions, indel_positions): output_indel_3 = [] output_indel_2 = ''.join(output_indel_2) records = output_indel_2.split('//') del records[0] ### Title manipulation titles = [] ls = records[0].split('\n') for i in range(6, len(ls), 1): a = ls[i].split(': ') titles.append(a[0]+': ') arr, i_details = [], [] for j in range(0, len(simple_indel_positions), 1): dat = [] if simple_indel_positions[j] == 's': lines = records[j].split('\n') a = lines[1].split(': ') no_of_indel = a[1] a = lines[2].split(': ') indel_location = a[1] a = lines[3].split(': ') b = a[1].split(' ') inde_size = b[0] del lines[:6] dat = [no_of_indel, indel_location, inde_size] i_details.append(dat) row = [] for line in lines: a = line.split(': ') seq = a[1].split(' ') if '-'*len(seq[1]) == seq[1]: row.append('0') elif not '-' in seq[1]: row.append('1') arr.append(row) arr = [[r[col] for r in arr] for col in range(len(arr[0]))] heading = '#NEXUS\nBEGIN DATA;\n DIMENSION NTAX=%s NCHAR=%s;\ \n FORMAT DATATYPE=SYMBOL "0 1";\ \n OPTIONS GAPMODE=MISSING;\n\nMATRIX' % (len(arr), len(arr[0])) output_indel_3.append(heading) r = getRuler(len(titles[0]), len(arr[0])).split('\n') ruler = '[' + r[0][1:] + ' '*(len(r[1])-len(r[0])) + ']\n[' + r[1][1:] + ']' output_indel_3.append(ruler) for k in range(0, len(titles), 1): line = titles[k] + str(''.join(arr[k])) output_indel_3.append(line) h1 = ';\nEND;\n\nBEGIN NOTES;\n[ Indel Number Alignment Position Indel Length ]\ \n[ ------------ ------------------ ------------ ]' output_indel_3.append(h1) for i_detail in i_details: line = '[ ' + i_detail[0] + ' '*(17-len(i_detail[0])) + i_detail[1] +\ ' '*(24-len(i_detail[1])) + i_detail[2] + ' '*(10-len(i_detail[2])) + ']' output_indel_3.append(line) output_indel_3.append('\nEND;') return output_indel_3 def getIndelCharacter(output_indel_2, simple_indel_positions, indel_positions): output_indel_3 = [] output_indel_2 = ''.join(output_indel_2) records = output_indel_2.split('//') del records[0] arr, i_details = [], [] for j in range(0, len(simple_indel_positions), 1): dat = [] if simple_indel_positions[j] == 's': lines = records[j].split('\n') a = lines[1].split(': ') no_of_indel = a[1] a = lines[2].split(': ') indel_location = a[1] a = lines[3].split(': ') b = a[1].split(' ') inde_size = b[0] del lines[:6] dat = [no_of_indel, indel_location, inde_size] i_details.append(dat) row = [] for line in lines: a = line.split(': ') seq = a[1].split(' ') if '-'*len(seq[1]) == seq[1]: row.append('0') elif not '-' in seq[1]: row.append('1') arr.append(row) return [[r[col] for r in arr] for col in range(len(arr[0]))] def indelExtraction(handle): pseudoalignments = genPseudoalignment(handle) conserved_block_profile = getConservedBlockProfile(pseudoalignments) indel_profile = getIndelProfile(conserved_block_profile) indel_positions = getIndelPositions(indel_profile) ruler = getRuler(len(handle[0][0]), len(handle[0][1])) ###################### ## OUTPUT Section ## ###################### if partial == 'True': output_indel_1 = genIndelAlignment(pseudoalignments, indel_profile, ruler, indel_positions) #output_indel_4 = genHomologAlignment(pseudoalignments, indel_profile, ruler, indel_positions) elif partial == 'False': output_indel_1 = genIndelAlignment(handle, indel_profile, ruler, indel_positions) #output_indel_4 = genHomologAlignment(handle, indel_profile, ruler, indel_positions) output_indel_4 = genHomologAlignment(handle, indel_profile, ruler, indel_positions) simple_indel_positions = search_for_simple_indels(pseudoalignments, inter_indels, indel_positions) output_indel_2 = genIndelLists(pseudoalignments, inter_indels, indel_positions, simple_indel_positions) if 's' in simple_indel_positions: output_indel_3 = getIndelMatrix(output_indel_2, simple_indel_positions, indel_positions) else: output_indel_3 = ['###################\n### NO MATRIX ###\n###################',\ 'SeqFIRE cannot generate indel matrix because there is no simple indel in the alignment!!'] filename = os.path.basename(infile).split('.') if output_mode == 1 or output_mode == 3: for out1 in output_indel_1: print (out1) print ('\n---SeqFIRE---\n') for out2 in output_indel_2: print (out2) print ('\n---SeqFIRE---\n') for out3 in output_indel_3: print (out3) print ('\n---SeqFIRE---\n') for out4 in output_indel_4: print (out4) if output_mode == 2 or output_mode == 3: ### Writing output1: Alignment with Indel Mask f1 = open(r'%s.txt' % (filename[0]), 'w') f1.write(output_indel_1[0]) f1.write(output_indel_1[1] + '\n') del output_indel_1[:2] for out1 in output_indel_1: f1.write(str(out1) + '\n') f1.close() ### Writing output2: Indel List f2 = open(r'%s.indel' % (filename[0]), 'w') f2.write(output_indel_2[0]) del output_indel_2[0] for out2 in output_indel_2: f2.write('\n' + str(out2)) f2.close() ### Writing output3: Indel Matrix f3 = open(r'%s_matrix.nex' % (filename[0]), 'w') for out3 in output_indel_3: f3.write('\n' + str(out3)) f3.close() ### Writing output4: Msked Alignment f4 = open(r'%s.nex' % (filename[0]), 'w') for out4 in output_indel_4: f4.write('\n' + str(out4)) f4.close() ######################################################################### ###### ###### ###### C O N S E R V E D B L O C K M O D U L E ###### ###### ###### ######################################################################### def getListOfSimilarityScores(handle): similarityList = [] for i in range(0, len(handle[0][1]), 1): aa_dict = {'A': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0,\ 'K': 0, 'L': 0, 'M': 0, 'N': 0, 'P': 0, 'Q': 0, 'R': 0, 'S': 0,\ 'T': 0, 'V': 0, 'W': 0, 'X': 0, 'Y': 0, '-': 0, '?': 0} for seq_list in handle: aa_dict[seq_list[1][i]] += 1 if p_matrix_2 == 'NONE': if (float(aa_dict['-']) * 100.0) / float(len(handle)) >= percent_accept_gap: similarityList.append('-') else: k = max(aa_dict, key=aa_dict.get) similarityList.append((float(aa_dict[k]) * 100.0) / float(len(handle))) else: pMatrix = pMatrices[p_matrix_2] if (float(aa_dict['-']) * 100.0) / float(len(handle)) >= percent_accept_gap: similarityList.append('-') else: c = [] for aa_sets in pMatrix: m = 0 for aa_set in aa_sets: m = m + int(aa_dict[aa_set]) c.append((float(m)*100.0)/float(len(handle))) similarityList.append(max(c)) return similarityList def getSimilarityProfile(similarity_list): similarityProfile = '' for i in similarity_list: if i == '-': similarityProfile = similarityProfile + i else: for onePair in percent_similarity: if onePair[0]<=i<=onePair[1]: similarityProfile = similarityProfile + 'H' break else: similarityProfile = similarityProfile + '.' return similarityProfile def getSimilarityBlocks(similarity_profile): similarityBlocks = similarity_profile for i in range(blocks, 0, -1): pattern = 'H' + '.'*i + 'H' replaceWith = 'H' + 'H'*i + 'H' while pattern in similarityBlocks: similarityBlocks = similarityBlocks.replace(pattern, replaceWith) for j in range(1, fuse+1, 1): pattern = '.' + 'H'*j + '.' replaceWith = '.' + '.'*j + '.' while pattern in similarityBlocks: similarityBlocks = similarityBlocks.replace(pattern, replaceWith) pattern = '.' + 'H'*j + '-' replaceWith = '.' + '.'*j + '-' while pattern in similarityBlocks: similarityBlocks = similarityBlocks.replace(pattern, replaceWith) pattern = '-' + 'H'*j + '.' replaceWith = '-' + '.'*j + '.' while pattern in similarityBlocks: similarityBlocks = similarityBlocks.replace(pattern, replaceWith) pattern = '-' + 'H'*j + '-' replaceWith = '-' + '.'*j + '-' while pattern in similarityBlocks: similarityBlocks = similarityBlocks.replace(pattern, replaceWith) return similarityBlocks def getInformationEntropy(handle): entropyList = [] for i in range(0, len(handle[0][1]), 1): aa_dict = {'A': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0,\ 'K': 0, 'L': 0, 'M': 0, 'N': 0, 'P': 0, 'Q': 0, 'R': 0, 'S': 0,\ 'T': 0, 'V': 0, 'W': 0, 'X': 0, 'Y': 0, '-': 0, '?': 0} for seq_list in handle: aa_dict[seq_list[1][i]] += 1 if (float(aa_dict['-']) * 100.0) / float(len(handle)) >= percent_accept_gap: entropyList.append('-') else: site_entropy = 0.0 for aa in amino_acid: proportion_of_any_amino_acid = float(aa_dict[aa])/float(len(handle)) if not proportion_of_any_amino_acid == 0.0: site_entropy = site_entropy + (proportion_of_any_amino_acid * log10(proportion_of_any_amino_acid)) site_entropy = site_entropy * (-1.0) entropyList.append(site_entropy) return entropyList def getMedian(numericValues): theValues = sorted(numericValues) if len(theValues) % 2 == 1: return theValues[(len(theValues)+1)//2-1] else: lower = theValues[len(theValues)//2-1] upper = theValues[len(theValues)//2] return (float(lower + upper)) / 2 def getEntropyProfile(entropy_list): datalists = [] for i in entropy_list: if not i == '-': datalists.append(i) median = getMedian(datalists) entropyProfile = '' for i in entropy_list: if i == '-': entropyProfile = entropyProfile + '-' elif i >= median*1.4826: entropyProfile = entropyProfile + '.' elif i < median*1.4826: entropyProfile = entropyProfile + 'E' #elif i >= median*1.4826 or i <= median*(-1.4826): entropyProfile = entropyProfile + '.' #elif i < median*1.4826 and i > median*(-1.4826): entropyProfile = entropyProfile + 'E' return entropyProfile def getEntropyBlocks(entropy_profile): entropyBlocks = entropy_profile for i in range(blocks, 0, -1): pattern = 'E' + '.'*i + 'E' replaceWith = 'E' + 'E'*i + 'E' while pattern in entropyBlocks: entropyBlocks = entropyBlocks.replace(pattern, replaceWith) for j in range(1, fuse+1, 1): pattern = '.' + 'E'*j + '.' replaceWith = '.' + '.'*j + '.' while pattern in entropyBlocks: entropyBlocks = entropyBlocks.replace(pattern, replaceWith) pattern = '.' + 'E'*j + '-' replaceWith = '.' + '.'*j + '-' while pattern in entropyBlocks: entropyBlocks = entropyBlocks.replace(pattern, replaceWith) pattern = '-' + 'E'*j + '.' replaceWith = '-' + '.'*j + '.' while pattern in entropyBlocks: entropyBlocks = entropyBlocks.replace(pattern, replaceWith) pattern = '-' + 'E'*j + '-' replaceWith = '-' + '.'*j + '-' while pattern in entropyBlocks: entropyBlocks = entropyBlocks.replace(pattern, replaceWith) return entropyBlocks def searchForConservedBlocks(similarity_blocks, entropy_blocks): conservedBlocks = '' if strick_combination == 'True': for i in range(0, len(similarity_blocks), 1): if similarity_blocks[i] == '-' and entropy_blocks[i] == '-': conservedBlocks = conservedBlocks + '-' elif similarity_blocks[i] == 'H' and entropy_blocks[i] == 'E': conservedBlocks = conservedBlocks + '#' else: conservedBlocks = conservedBlocks + '.' elif strick_combination == 'False': for i in range(0, len(similarity_blocks), 1): if similarity_blocks[i] == '-' and entropy_blocks[i] == '-': conservedBlocks = conservedBlocks + '-' elif similarity_blocks[i] == 'H' or entropy_blocks[i] == 'E': conservedBlocks = conservedBlocks + '#' else: conservedBlocks = conservedBlocks + '.' return conservedBlocks def genConservedAlignment(handle, conserved_blocks, ruler, indel_profile): output_conserved_1 = [] heading = '%s\n# OUTPUT: PROTEIN ALIGNMENT WITH CONSERVED BLOCK PROFILE\n# NOTE: # in the reference line means conserved block positions.\n' % (prog_title) output_conserved_1.append(heading) output_conserved_1.append(ruler) for i in handle: line = i[0] + i[1] output_conserved_1.append(line) line = 'Conserved Blocks' + ' ' * (len(handle[0][0])-16) + conserved_blocks output_conserved_1.append(line) if indel_profile == '': pass else: line = 'Indel Profile ' + ' ' * (len(handle[0][0])-16) + indel_profile output_conserved_1.append(line) return output_conserved_1 def genFastaWithConservedProfile(handle, conserved_blocks): output_conserved_2 = [] for record in handle: output_conserved_2.append('>' + record[0][:-2]) output_conserved_2.append(record[1]) output_conserved_2.append('>Conserved blocks') output_conserved_2.append(conserved_blocks) return output_conserved_2 def genFastaWithConservedBlockOnly(handle, conserved_blocks): output_conserved_3, conserved_index = [], [] # Find the conserved block positions ref = '-' + conserved_blocks + '-' for i in range(0, len(ref)-1, 1): current_position = ref[i] next_position = ref[i+1] if not current_position == '#' and next_position == '#': conserved_index.append(i) elif current_position == '#' and not next_position == '#': conserved_index.append(i) for record in handle: output_conserved_3.append('>' + record[0][:-2]) seq = '' for i in range(0, len(conserved_index), 2): seq = seq + record[1][conserved_index[i]:conserved_index[i+1]] output_conserved_3.append(seq) return output_conserved_3 def genNexusWithConservedProfile(handle, conserved_blocks, simple_indel_positions, indel_positions, indel_profile): output_conserved_4 = [] conserved_index = [] charset = ' charset conRes =' # Prepare conRes ref = '-' + conserved_blocks + '-' for i in range(0, len(ref)-1, 1): current_position = ref[i] next_position = ref[i+1] if not current_position == '#' and next_position == '#': conserved_index.append(i) elif current_position == '#' and not next_position == '#': conserved_index.append(i) for i in range(0, len(conserved_index), 2): charset = charset + ' ' + str(conserved_index[i]+1) + '-' + str(conserved_index[i+1]) head = '''#NEXUS BEGIN DATA; \tDIMENSIONS NTAX=%s NCHAR=%s; \tFORMAT MISSING=? DATATYPE=PROTEIN GAP=- EQUATE="0=K 1=D"; \tOPTIONS GAPMODE=MISSING; MATRIX''' % (len(handle),len(handle[0][1])) output_conserved_4.append(head) r = getRuler(len(handle[0][0]), len(handle[0][1])).split('\n') ruler = '[ ' + r[0][1:] + ' '*(len(r[1])-len(r[0])) + ']\n[ ' + r[1][1:] + ']' output_conserved_4.append(ruler) for i in handle: line = i[0] + ' ' + i[1] output_conserved_4.append(line) line = '[Conserved Block Profile' + ' ' * (len(handle[0][0])-23) + conserved_blocks + ']' output_conserved_4.append(line) indelChar = '' for j in range(0, len(indel_positions), 2): indelChar = indelChar + ' ' + str(indel_positions[j]) + '-' + str(indel_positions[j+1]) if combine_with_indel == 'True': line = '[Indel profile' + ' ' * (len(handle[0][0])-13) + indel_profile + ']' output_conserved_4.append(line) line = ';\nEND;\n\nBEGIN SETS;' output_conserved_4.append(line) charset = charset + ';\n charset indels =' + indelChar charset = charset + ';\n charset indelMatrix = ' + str(len(handle[0][1])-simple_indel_positions.count('s')+1) + '-' + str(len(handle[0][1])) + ';\nEND;' output_conserved_4.append(charset) #line = '\nBEGIN CHARACTER;\n eliminate ' #output_conserved_4.append(line) h1 = 'END;\n\nBEGIN NOTES;\n[ Indel Number Alignment Position Indel Length ]\n[ ------------ ------------------ ------------ ]' output_conserved_4.append(h1) for i, status in enumerate(simple_indel_positions): if status == 's': if indel_positions[0]+1 == indel_positions[1]: align_pos = str(indel_positions[0]+1) else: align_pos = str(indel_positions[0]+1) + '-' + str(indel_positions[1]) indel_length = indel_positions[1] - indel_positions[0] line = '[ ' + str(i+1) + ' '*(17-len(str(i+1))) + align_pos + ' '*(24-len(align_pos)) + str(indel_length) + ' '*(9-len(str(indel_length))) + ']' output_conserved_4.append(line) del indel_positions[0:2] elif status == 'c': del indel_positions[0:2] output_conserved_4.append('\nEND;') elif combine_with_indel == 'False': line = ';\nEND;\n\nBEGIN SETS;' output_conserved_4.append(line) charset = charset + ';\n charset indels =' + indelChar + ';\nEND;' output_conserved_4.append(charset) #line = '\n\nBEGIN CHARACTER;\n eliminate ' #output_conserved_4.append(line) #line = 'END;' #output_conserved_4.append(line) return output_conserved_4 def genNexusWithConservedBlockOnly(output_conserved_3): output_conserved_5 = [] head = '#NEXUS\nBEGIN DATA;\n DIMENSIONS NTAX=%s NCHAR=%s;\ \n FORMAT MISSING=? DATATYPE=PROTEIN GAP=- EQUATE="0=K 1=D";\ \n OPTIONS GAPMODE=MISSING;\n\nMATRIX' % (len(output_conserved_3)//2, len(output_conserved_3[1])) output_conserved_5.append(head) r = getRuler(len(output_conserved_3[0]), len(output_conserved_3[1])).split('\n') ruler = '[ ' + r[0][1:] + ' '*(len(r[1])-len(r[0])) + ']\n[ ' + r[1][1:] + ']' output_conserved_5.append(ruler) for i in range(0, len(output_conserved_3), 2): line = output_conserved_3[i][1:] + ': ' + output_conserved_3[i+1] output_conserved_5.append(line) line = ';\nEND;' output_conserved_5.append(line) return output_conserved_5 def conservedBlockExtraction(handle): similarity_list = getListOfSimilarityScores(handle) similarity_profile = getSimilarityProfile(similarity_list) similarity_blocks = getSimilarityBlocks(similarity_profile) entropy_list = getInformationEntropy(handle) entropy_profile = getEntropyProfile(entropy_list) entropy_blocks = getEntropyBlocks(entropy_profile) conserved_blocks = searchForConservedBlocks(similarity_blocks, entropy_blocks) ############################################ ## Extraction of Conserved Blocks ## ############################################ if combine_with_indel == 'True': pseudoalignments = genPseudoalignment(handle) conserved_block_profile = getConservedBlockProfile(pseudoalignments) indel_profile = getIndelProfile(conserved_block_profile) indel_positions = getIndelPositions(indel_profile) simple_indel_positions = search_for_simple_indels(pseudoalignments, inter_indels, indel_positions) output_indel_2 = genIndelLists(pseudoalignments, inter_indels, indel_positions, simple_indel_positions) if 's' in simple_indel_positions: indel_in_matrix = getIndelCharacter(output_indel_2, simple_indel_positions, indel_positions) for i in range(0, len(indel_in_matrix), 1): handle[i][1] = handle[i][1] + ''.join(indel_in_matrix[i]) ruler = getRuler(len(handle[0][0]), len(handle[0][1])) output_conserved_1 = genConservedAlignment(handle, conserved_blocks, ruler, indel_profile) output_conserved_2 = genFastaWithConservedProfile(handle, conserved_blocks) output_conserved_3 = genFastaWithConservedBlockOnly(handle, conserved_blocks) output_conserved_4 = genNexusWithConservedProfile(handle, conserved_blocks, simple_indel_positions, indel_positions, indel_profile) output_conserved_5 = genNexusWithConservedBlockOnly(output_conserved_3) elif combine_with_indel == 'False': pseudoalignments = genPseudoalignment(handle) conserved_block_profile = getConservedBlockProfile(pseudoalignments) indel_profile = getIndelProfile(conserved_block_profile) indel_positions = getIndelPositions(indel_profile) ruler = getRuler(len(handle[0][0]), len(handle[0][1])) output_conserved_1 = genConservedAlignment(handle, conserved_blocks, ruler, '') output_conserved_2 = genFastaWithConservedProfile(handle, conserved_blocks) output_conserved_3 = genFastaWithConservedBlockOnly(handle, conserved_blocks) output_conserved_4 = genNexusWithConservedProfile(handle, conserved_blocks, '', indel_positions, '') output_conserved_5 = genNexusWithConservedBlockOnly(output_conserved_3) ###################### ## OUTPUT Section ## ###################### filename = os.path.basename(infile).split('.') if output_mode == 1 or output_mode == 3: for out1 in output_conserved_1: print (out1) print ('\n---SeqFIRE---\n') for out2 in output_conserved_2: print (out2) print ('\n---SeqFIRE---\n') for out3 in output_conserved_3: print (out3) print ('\n---SeqFIRE---\n') for out4 in output_conserved_4: print (out4) print ('\n---SeqFIRE---\n') for out5 in output_conserved_5: print (out5) if output_mode == 2 or output_mode == 3: ### Writing output1: Alignment with Indel Mask f1 = open(r'%s.txt' % (filename[0]), 'w') f1.write(output_conserved_1[0]) f1.write(output_conserved_1[1] + '\n') del output_conserved_1[:2] for out1 in output_conserved_1: f1.write(str(out1) + '\n') f1.close() ### Writing output2: Alignment in Fasta format + Conserved block profile f2 = open(r'%s.fasta' % (filename[0]), 'w') for out2 in output_conserved_2: f2.write('\n' + str(out2)) f2.close() ### Writing output3: Alignment in Fasta format f3 = open(r'%s_2_short.fasta' % (filename[0]), 'w') for out3 in output_conserved_3: f3.write('\n' + str(out3)) f3.close() ### Writing output4: Alignment in NESXUS format + Conserved block profile f4 = open(r'%s_2.nex' % (filename[0]), 'w') for out4 in output_conserved_4: f4.write('\n' + str(out4)) f4.close() ### Writing output4: Alignment in NESXUS format + Conserved block profile f5 = open(r'%s_2_short.nex' % (filename[0]), 'w') for out5 in output_conserved_5: f5.write('\n' + str(out5)) f5.close() ######################################################################### ###### ###### ###### S E Q U E N C E V E R I F I C A T I O N ###### ###### ###### ######################################################################### ######################## # RUN BEFORE ANALYZING # ######################## #This function check whether each input sequence is in FASTA format def checkSeqFormat(inputText): if '>' not in inputText: return {'Fatal':True, 'error_messages':['The sequence is not in FASTA format']} splitInput = inputText.split('>') if splitInput[0]!='': return {'Fatal':True, 'error_messages':['The sequence is not in FASTA format']} del splitInput[0] #print(splitInput) for oneInput in splitInput: if oneInput!=oneInput.strip(' '): return {'Fatal':True, 'error_messages':['The sequence is not in FASTA format']} return True #This function check whether the input has been prepared by seqFIREprep def checkPrepped(inputText): if '==seq==' not in inputText: return {'Fatal':True, 'error_messages':['The sequence is not properly prepped']} splitInput = inputText.split('==seq==') if splitInput[0]!='': return {'Fatal':True, 'error_messages':['The sequence is not properly prepped']} del splitInput[0] for oneInput in splitInput: if checkSeqFormat(oneInput)!=True: return {'Fatal':True, 'error_messages':['One of the sequence is not in FASTA format']} return True ####################### # RUN AFTER ANALYZING # ####################### ''' Input for these functions is "seqList" [[,], [,], ...] ''' #This function check whether each input sequence is DNA or amino acid sequence ''' Input: seqList and user selected sequence (selectedType) Output: True if no error occurred, error list otherwise A sequence is not both DNA and amino acid sequence when 1. The sequence consists of "J", "O", "U", "Z" A sequence predicted as either DNA or amino acid sequence otherwise Some letter in the sequence is extracted into 2 groups 1. notDNA group: "E", "F", "I", "L", "P", "Q", "X" 2. notProtein group: "B" The sequence is then classified with this criteria: 1. If the sequence consists of both notDNA and notProtein group -> Cannot determine 2. If the sequence NOT consists of both notDNA and notProtein group If at least 90% of the sequence is "A", "T", "C", "G" -> DNA Otherwise -> Protein 3. If the sequence only consists of notDNA group -> Protein 4. If the sequence only consists of notProtein group -> DNA The predicted sequence type (Verdict) is checked with user's selected sequence type Error will be raised when the verdict and selected sequence type does not match ''' def checkSeqType(seqList,selectedType): verdict = "" errormsg = [] for oneSeq in seqList: letter_dict = {'A': 0, 'B': 0, 'C': 0, 'D': 0, 'E': 0, 'F': 0, 'G': 0, 'H': 0, 'I': 0,'K': 0, 'L': 0, 'M': 0, 'N': 0, 'P': 0, 'Q': 0, 'R': 0, 'S': 0, 'T': 0, 'V': 0, 'W': 0, 'X': 0, 'Y': 0, '-': 0, '?': 0} #Exclude J, O, U, Z for oneLetter in oneSeq[1]: try: letter_dict[oneLetter] = letter_dict[oneLetter] + 1 except: errormsg.append("FATAL ERROR: Detected invalid character in sequence %s" % (oneSeq[0])) break #Making prediction whether the sequence is DNA or RNA notDNAScore = letter_dict['E']+letter_dict['F']+letter_dict['I']+letter_dict['L']+letter_dict['P']+letter_dict['Q']+letter_dict['X'] notProteinScore = letter_dict['B'] if notDNAScore==0 and notProteinScore==0: DNAScore = letter_dict['A']+letter_dict['T']+letter_dict['C']+letter_dict['G'] if DNAScore>=0.9*(sum(letter_dict.values())-letter_dict['-']-letter_dict['?']): verdict='DNA' else: verdict='Protein' elif notDNAScore>0 and notProteinScore>0: errormsg.append("FATAL ERROR: Detected invalid character in sequence %s" % (oneSeq[0])) continue elif notDNAScore>0: verdict='Protein' elif notProteinScore>0: verdict='DNA' #Test verdict with selected sequence type if verdict=="DNA" and selectedType == "Protein": errormsg.append("The sequence %s is detected as DNA" % (oneSeq[0])) elif verdict=="Protein" and selectedType =="DNA": errormsg.append("The sequence %s is detected as amino acid" % (oneSeq[0])) return errormsg if errormsg else True #This function check if the input sequences are multiple aligned ''' Input: seqList Output: True if seqList is MSA, False otherwise A sequence is MSA when the length of seqList is more than or equal to 2 All sequences must have the same length, or it will not be a valid MSA ''' def checkMultipleSeq(seqList): if len(seqList)>=2: seqLength = len(seqList[0][1]) for oneSeq in seqList: if len(oneSeq[1])!=seqLength: return ["FATAL ERROR: Inequal sequence length detected"] return True else: return ["FATAL ERROR: Single alignment detected"] #This function check the alignment quality ''' Input: seqList Output: True if the quality of MSA is valid, error list otherwise Validity of a sequence is evaluated from 2 criteria: Row validity and Column validity 1. Row validity - A sequence with at least 40% CONTINUOUS GAP at head or tail will be considered invalid ''' def checkMSAQuality(seqList): errormsg = [] for oneSeq in seqList: headGap = 0 tailGap = 0 if oneSeq[1][0]=='-': for oneLetter in oneSeq[1]: if oneLetter=='-': headGap = headGap+1 else: break if oneSeq[1][-1]=='-': for oneLetter in oneSeq[1][::-1]: if oneLetter=='-': tailGap = tailGap+1 else: break biggestGap = headGap if headGap>tailGap else tailGap if biggestGap>=0.4*len(oneSeq[1]): errormsg.append('The sequence %s has at least' % oneSeq[0] + ' continuous gap at head or tail') return errormsg if errormsg else True #This function combines all 4 checking functions into one ''' Input: seqList Output: True if sequences are ready to analyze ''' def checkReadiness(seqList, selectedType): errors = [] errorDict = {'Fatal':False,'error_messages':[]} seqTypeResult = checkSeqType(seqList, selectedType) multipleSeqResult = checkMultipleSeq(seqList) MSAQualityResult = checkMSAQuality(seqList) if seqTypeResult!=True: errors.extend(seqTypeResult) if multipleSeqResult!=True: errors.extend(multipleSeqResult) if MSAQualityResult!=True: errors.extend(MSAQualityResult) for oneError in errors: if oneError.startswith('FATAL ERROR: '): errorDict['Fatal']=True errorDict['error_messages'].append(oneError) if errorDict['error_messages']!=[]: return errorDict else: return True def showError(errorDict): print('***************\n* ERROR FOUND *\n***************') if errorDict['Fatal']==True: raise Exception(errorDict['error_messages'][0]) else: for oneError in errorDict['error_messages']: print(oneError) print("\n") ######################################################################### ###### ###### ###### S e q F I R E ' s M A I N P R O G R A M ###### ###### ###### ######################################################################### if sys.argv[1:]==[]: usage() sys.exit(2) try: opts, args = getopt.getopt(sys.argv[1:], "h:i:a:c:d:j:g:k:b:t:p:s:f:r:e:m:o:") except (getopt.GetoptError): print("Invalid arguments") #usage() sys.exit(2) for opt, arg in opts: if opt == "-h": usage() sys.exit() if opt == "-i": infile = arg if opt == "-a": analysis_mode = int(arg) if opt == "-c": similarity_threshold = float(arg) if opt == "-d": percent_similarity = float(arg) if opt == "-j": percent_accept_gap = float(arg) if opt == "-g": p_matrix = arg if opt == "-k": p_matrix_2 = arg if opt == "-b": inter_indels = int(arg) if opt == "-s": blocks = int(arg) if opt == "-f": fuse = int(arg) if opt == "-r": strick_combination = str(arg) if opt == "-e": combine_with_indel = str(arg) if opt == "-m": multidata = int(arg) if opt == "-o": output_mode = int(arg) if opt == '-t': seqType = arg ################################### ## A N A L Y S I S Z O N E ## ################################### if multidata == 1: f = open(r'%s' % (infile), 'r') record = f.read() f.close() checkSeqResult = checkSeqFormat(record) if checkSeqResult == True: handle = parseFasta(record) else: showError(checkSeqResult) readinessResult = checkReadiness(handle,seqType) if readinessResult != True: showError(readinessResult) if analysis_mode == 1: indelExtraction(handle) ### INDEL REGION MODULE ### elif analysis_mode == 2: conservedBlockExtraction(handle) ### CONSERVED BLOCK MODULE ### elif multidata == 2: f = open(r'%s' % (infile), 'r') records = f.read() f.close() checkPreppedResult = checkPrepped(records) if checkPreppedResult == True: records = records.split('==seq==') del records[0] else: showError(checkPreppedResult) for record in records: a = record.split('==fire==') handle = parseFasta(a[1]) readinessResult = checkReadiness(handle,seqType) if readinessResult != True: showError(readinessResult) if output_mode == 1 or output_mode == 3: print ('==seq==%s==fire==' % infile) elif output_mode == 2 or output_mode == 3: f = open(r'outfile.txt', 'a') f.write('\n==seq==%s==fire==\n' % infile) f.close() if analysis_mode == 1: indelExtraction(handle) ### INDEL REGION MODULE ### elif analysis_mode == 2: conservedBlockExtraction(handle) ### CONSERVED BLOCK MODULE ###