# Bioinformatics

## DNA Alignment

This takes two DNA sequences and produces the optimal alignment. This is done by filling out a backtracking matrix.

This has a simple piece of driver code which loops through all the cells of the backtrack matrix

``````for a in range(1,len(seq1)+1):
for b in range(1,len(seq2)+1):
fill_cell(m,backtrack,seq1,seq2,a,b)
``````

The `fill_cell` function then determines the entry in the backtracking matrix and the maximum score of a matching

``````def fill_cell(m,backtrack,seq1,seq2,a,b):
# Here m is the matrix and [b,a] is the location of the cell to fill out
# Diagonal
max=score(a,b,seq1,seq2)+m[b-1,a-1]
pos='D'
# Up
temp=m[b-1,a]-2
if temp>max:
max=temp
pos='U'
# Left
temp=m[b,a-1]-2
if temp>max:
max=temp
pos='L'
m[b,a]=max
backtrack[b,a]=pos
return
``````

The rules for scoring were given in the assignment, and are implemented using the `score` function

``````def score(a,b,seq1,seq2):
# This function finds the score of the matching
index_a=a-1
index_b=b-1
if seq1[index_a]==seq2[index_b]:
if seq1[index_a]=='A':
return 4
if seq1[index_a]=='C':
return 3
if seq1[index_a]=='G':
return 2
if seq1[index_a]=='T':
return 1
else:
print(str(seq1[index_a])+str(seq2[index_b])+str(index_a)+str(index_b))
else:
return -3
``````

The best alignment is then generated using the following function by moving through the backtracking matrix

``````def gen_seq(backtrack,str1,str2):
coord=(len(str2),len(str1))
matchstring1=''
matchstring2=''
while coord!=(0,0):
if backtrack[coord]=='D':
coord=(coord[0]-1,coord[1]-1)
matchstring1=matchstring1+str1[-1]
str1=str1[:-1]
matchstring2 = matchstring2 + str2[-1]
str2 = str2[:-1]
if backtrack[coord]=='U':
coord=(coord[0]-1,coord[1])
matchstring1 = matchstring1 + '-'
matchstring2 = matchstring2 + str2[-1]
str2 = str2[:-1]
if backtrack[coord]=='L':
coord=(coord[0],coord[1]-1)
matchstring1 = matchstring1 + str1[-1]
str1 = str1[:-1]
matchstring2 = matchstring2 + '-'
matchstring1=matchstring1[::-1]
matchstring2 = matchstring2[::-1]
return [matchstring1,matchstring2]
``````

### Drawing a phylogenetic tree

This task is to generate a phylogenetic tree based on an input distance matrix

For brevity I have removed the code which formats the input from the file so that it can be worked on along with the code to output the tree.

``````def WPGMA(filename):
while table.shape != (2, 2):
# Find the minimum value
minval = np.min(table[np.nonzero(table)])
# Find its coordinates
itemindex = np.where(table == minval)
# Merge the two species corresponding to those coordinates
table = mergespecies(table, itemindex[0][0], itemindex[0][1], names, names2)
nametable=[str(name) for name in names2]
stack1=np.array(nametable)
stacktotal=np.vstack((stack1,table))
print(stacktotal)
``````

And the function to merge species is as follows:

``````def mergespecies(table, a, b, names, names2):
depth = lambda L: (isinstance(L, list) and (max(map(depth, L)) + 1) if L else 1) or 0
if depth(names[a])<depth(names[b]):
names[a],names[b]=names[b],names[a]
# print("Merge:" + str(names[a]) + " and " + str(names[b]))
graphmerge(names[a],names[b])
names2[a]=names2[a]+names2[b]
names2.remove(names2[b])
sublist = [names[a], names[b]]
names.remove(names[b])
names.remove(names[a])
names.insert(a, sublist)
column1 = table[:, a:a + 1]
column2 = table[:, b:b + 1]
# Combine the two columns
combine = np.hstack((column1, column2))
# Get the mean of all the rows
combine = combine.reshape(-1, 2).mean(axis=1).reshape(combine.shape[0], -1)
# Delete the rows corresponding to the two species
combine = np.delete(combine, [a, b], 0)
# Delete the rows/columns corresponding to the two species from the main table
table = np.delete(table, [a, b], 0)
table = np.delete(table, [a, b], 1)
# Append the amended column to the main table
combine = np.transpose(combine)
table = np.insert(table, min(a, b), combine, axis=1)
# Add a zero the the column and transpose it
combine = np.insert(combine, min(a, b), [0], axis=1)
# Add this row to the bottom of the main table
table = np.insert(table, min(a, b), combine, axis=0)
# Return the table back to the main function
return table
``````