【ROSALIND】【练Python，学生信】41 构建距离矩阵

如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面  谢谢配合~

题目：

构建距离矩阵（Creating a Distance Matrix）

Given: A collection of n (n≤10) DNA strings s1,…,sn of equal length (at most 1 kbp). Strings are given in FASTA format.

所给：不多于十条DNA序列，长度都相等（最多1kb），以FASTA格式给出。

Return: The matrix D corresponding to the p-distance dp on the given strings. As always, note that your answer is allowed an absolute error of 0.001.

需得：所给序列的距离矩阵，允许绝对误差不超过0.001。

测试数据

>Rosalind_9499

TTTCCATTTA

>Rosalind_0942

GATTCATTTC

>Rosalind_6568

TTTCCATTTT

>Rosalind_1833

GTTCCATTTA

测试输出

0.00000 0.40000 0.10000 0.10000

0.40000 0.00000 0.40000 0.30000

0.10000 0.40000 0.00000 0.20000

0.10000 0.30000 0.20000 0.00000

生物学背景

        构造进化树有多种方法，各有优缺点，其中一种是基于距离矩阵的方法。即，用不同类之间的距离来构建树。计算距离的方法也有多种，这里采用Hamming距离，有关Hamming距离的介绍请见06 DNA序列Hamming距离的计算。

        对两个长度相等的DNA序列s1和s2，计算遗传距离最简单的方法是算p距离（p-distance），记为dp(s1,s2)，是s1和s2之间不同的字符数占字符总数的比例。对n个长度相等的DNA序列s1，s2…sn，其距离矩阵为Di,j=d(si,sj)。

思路

        这道题我直接用三层循环实现，前两层将序列两两取出，在最后一层循环依次比较。

代码

def readfasta(lines):
    """读入fasta格式文件的函数"""
    seq = []
    index = []
    seqplast = ""
    numlines = 0
    for i in lines:
        if '>' in i:
            index.append(i.replace("\n", "").replace(">", ""))
            seq.append(seqplast.replace("\n", ""))
            seqplast = ""
            numlines += 1
        else:
            seqplast = seqplast + i.replace("\n", "")
            numlines += 1
        if numlines == len(lines):
            seq.append(seqplast.replace("\n", ""))
    seq = seq[1:]
    return index, seq


f = open('input.txt', 'r')
lines = f.readlines()
f.close()

[index, seq] = readfasta(lines)
length = len(seq[0]) # 每个序列的字符数

result = []
for i in range(0, len(seq)): # 以下两个循环实现序列两两取出
    tep = []
    for j in range(0, len(seq)):
        value = 0
        for k in range(length): # 最内层循环实现两序列间字符依次比较
            if seq[i][k] != seq[j][k]:
                value += 1
        tep.append(value/length)
    result.append(tep)

f = open('output.txt', 'a')
for i in range(len(result)): # 按所给格式输出
    for j in range(len(result[i])):
        f.write(str(result[i][j]))
        f.write(' ')
    f.write('\n')

f.close()