Write a function printBRH$($geneName$,$allScoresD$)$

This function takes a gene name, and a dictionary of memoAlignScore scores as input. It finds and prints the best reciprocal hit for geneName. If there is no best reciprocal hit it returns without printing anything.

def memoAlignScore$($S$1,$ S$2,$ gap, substitutionMatrix, memo$)$:

$\text{'}\text{'}\text{'}$computes the alignment score using a a memoized version of alignScore,

memoAlignScore$($S$1,$ S$2,$ gap, substitutionMatrix,memo$),$ where memo is a

dictionary.$\text{'}\text{'}\text{'}$

# Check if result is already in memo

if $($S$1,$ S$2)$ in memo:

return memo$[($S$1,$ S$2)]$

# Base cases

if S$1==""$: return gap $*$ len$($S$2)$

elif S$2==""$: return gap $*$ len$($S$1)$

# Recursive cases

option$1=$ substitutionMatrix$[($S$1[0],$ S$2[0])]+$ memoAlignScore$($S$1[1$:$],$ S$2[1$:$],$ gap, substitutionMatrix, memo$)$

option$2=$ gap $+$ memoAlignScore$($S$1[1$:$],$ S$2,$ gap, substitutionMatrix, memo$)$

option$3=$ gap $+$ memoAlignScore$($S$1,$ S$2[1$:$],$ gap, substitutionMatrix, memo$)$

# Compute max score

best$\_$score $=$ max$($option$1,$ option$2,$ option$3)$

# Store result in memo

memo$[($S$1,$ S$2)]=$ best$\_$score

return best$\_$score

A simple way to find the best reciprocal hit for geneName is to call closestMatch to get its best matching gene in the other species. Then call closestMatch on this best match from the other species. If the output after the round trip matches geneName, then we've found a best reciprocal hit.

def closestMatch$($geneName$,$ allScoresD$)$:

$\text{'}\text{'}\text{'}$Given a gene name and a dictionary of alignment scores, closestMatch

returns the protein from the other species which is most similar

$($has the highest alignment score$).\text{'}\text{'}\text{'}$

# Initialize variables to keep track of the highest score and best matching protein

max$\_$score $=-$float$(\text{'}$inf$\text{'})$ # Use negative infinity to ensure any score will be higher

best$\_$match $=$ None

# Loop through each key$-$value pair in the dictionary

for $($geneA$,$ geneB$),$ score in allScoresD.items$()$:

# Check if geneName is in the current tuple $($geneA$,$ geneB$)$

if geneName in $($geneA$,$ geneB$)$:

# Determine the other gene in the pair

other$\_$gene $=$ geneB if geneA $==$ geneName else geneA

# Update max$\_$score and best$\_$match if this score is higher

if score $>$ max$\_$score:

max$\_$score $=$ score

best$\_$match $=$ other$\_$gene

return best$\_$match

For each best reciprocal hit we will print chromosome, start position, and gene name for the gene in the first species, and then for the gene in the second species. For example:

$>>>$ allScoresD$=$allScores$($sampleHumanGeneList$,$sampleChickenGeneList$)$

$>>>$ printBRH$(\text{'}$c$8\text{'},$allScoresD$)$

chr$243123243$ c$8---$ chr$345016733$ h$17$

$>>>$ printBRH$(\text{'}$h$7\text{'},$allScoresD$)$

$>>>$ printBRH$(\text{'}$h$17\text{'},$allScoresD$)$

chr$345016733$ h$17---$ chr$243123243$ c$8$