Consensus and Profile

Finding a Most Likely Common Ancestor

In “Counting Point Mutations”, we calculated the minimum number of symbol mismatches between two strings of equal length to model the problem of finding the minimum number of point mutations occurring on the evolutionary path between two homologous strands of DNA. If we instead have several homologous strands that we wish to analyze simultaneously, then the natural problem is to find an average-case strand to represent the most likely common ancestor of the given strands.

Problem

A matrix is a rectangular table of values divided into rows and columns. An m×n matrix has m rows and n columns. Given a matrix A, we write Ai,j to indicate the value found at the intersection of row i and column j.

Say that we have a collection of DNA strings, all having the same length n. Their profile matrix is a 4 × n matrix P in which P1,j represents the number of times that ‘A’ occurs in the jth position of one of the strings, P2,j represents the number of times that C occurs in the j-th position, and so on.

A consensus string c is a string of length n formed from our collection by taking the most common symbol at each position; the j-th symbol of c therefore corresponds to the symbol having the maximum value in the j-th column of the profile matrix. Of course, there may be more than one most common symbol, leading to multiple possible consensus strings.

Given

A collection of at most 10 DNA strings of equal length (at most 1 kbp) in FASTA format.

Return

A consensus string and profile matrix for the collection (If several possible consensus strings exist, then you may return any one of them).

Solution

Let’s define a function that would find out the locations of substrings:

def read_fasta(fp):
    name, seq = None, []
    for line in fp:
        line = line.rstrip()
        if line.startswith(">"):
            if name: yield (name, ''.join(seq))
            name, seq = line, []
        else:
            seq.append(line)
    if name: yield (name, ''.join(seq))

Then, read the dataset and run the function:

data_list = []
file = open('main.out','w')
with open('rosalind_cons.txt') as fp: #replace filename with yours
    for name, seq in read_fasta(fp):
        data_list.append(seq)
length = len(data_list)

L = len(seq)

P = [[0 for x in xrange(L)] for y in xrange(4)] 

Q = ['A','C','G','T']

for x in range(L):
    for y in range(4):
        for z in range(length):
            P[y][x] = P[y][x] + data_list[z][x].count(Q[y])

domi = ""
for x in range(L):
    MAX = 0
    for y in range(4):  
        if P[y][x]>=P[MAX][x]:
            MAX = y       
            
    if MAX == 0:
        domi = domi+'A'
    elif MAX ==1:
        domi = domi+'C'
    elif MAX ==2:
        domi = domi+'G'
    elif MAX ==3:
        domi = domi+'T'

And finally, print the output:

file.write('%s\n%s\n%s\n%s\n%s' %(
    domi,'A: '+str(P[0]).strip('[]').replace(',',''),'C: '
    +str(P[1]).strip('[]').replace(',',''),'G: '
    +str(P[2]).strip('[]').replace(',',''),'T: '
    +str(P[3]).strip('[]').replace(',','')))