生信刷题之ROSALIND——Part 5 (PERM, PRTM, REVP)

公众号搜索《生信er》，内容更多，更精彩~

1、Enumerating Gene Orders

Problem

点突变可以在同一物种的生物群体中产生变化，但它们缺乏创造和分化整个物种的能力。基因组重排（Genome rearrangement）是一种主要的基因组突变，通常是由减数分裂或有丝分裂后细胞分裂的错误引起的。这些对染色体结构的大规模改变几乎总是有害的，通常会导致发育中的生物体死亡或不育，但在极少数情况下，它们提供了显著的优势。

因为影响物种进化的重排很少发生，两个密切相关的物种将有非常相似的基因组。因此，为了简化两个这样的基因组的比较，研究人员首先从物种中找出相似的DNA片段，称为synteny blocks。随着时间的推移，基因组重排创造了这些块，并将它们在两个基因组之间来回移动（通常将这些块分离到不同的染色体上）。

因此，我们可以用一个正整数来标记每个synteny blocks；当比较两个物种的基因组/染色体时，我们只需要指定其编号的synteny blocks的顺序。

Given： 一个正整数 n≤7。

Return： 长度为n的排列的总数，后面是所有排列组合（以任意顺序排列）。

Sample Dataset

3

Sample Output

6
1 2 3
1 3 2
2 1 3
2 3 1
3 1 2
3 2 1

example

类似于数据结构中的全排列
解决方法：递归法
以{1,2,3}为例，它的排列是：
以1开头，后面接着{2,3}的全排列，
以2开头，后面接着{1,3}的全排列，
以3开头，后面接着{1,2}的全排列。

Code

# Enumerating Gene Orders
from itertools import permutations


# permutations 返回可迭代对象的所有数学或者字符的全排列方式
def get_permutations(n):
    seq_list = []
    for i in range(n):
        seq_list.append(i + 1)
    seqs = list(permutations(seq_list, n))
    print(len(seqs))
    for seq in seqs:
        print(" ".join(str(m) for m in seq))


get_permutations(3)

with open("rosalind_perm.txt", "r") as f:
    num = f.read().strip()
    num = int(num)
    get_permutations(num)

Output

6
1 2 3
1 3 2
2 1 3
2 3 1
3 1 2
3 2 1

720
1 2 3 4 5 6
1 2 3 4 6 5
...........
6 5 4 1 3 2
6 5 4 2 1 3
6 5 4 2 3 1
6 5 4 3 1 2
6 5 4 3 2 1

2、Calculating Protein Mass

Problem

在加权字母表中，每个符号都被赋予一个正实数，称为权重。由加权字母表中符号形成的字符串称为加权字符串，其权重等于其符号的权重之和。
分配给20个符号的氨基酸字母的每个成员的标准重量是相应的氨基酸的单异构质量。

Given： 一个蛋白质串P，长度最多为1000 aa。

Return： P的总重量。

# Monoisotopic mass table
A   71.03711
C   103.00919
D   115.02694
E   129.04259
F   147.06841
G   57.02146
H   137.05891
I   113.08406
K   128.09496
L   113.08406
M   131.04049
N   114.04293
P   97.05276
Q   128.05858
R   156.10111
S   87.03203
T   101.04768
V   99.06841
W   186.07931
Y   163.06333 

Sample Dataset

SKADYEK

Sample Output

821.392

Code

# Monoisotopic mass table

def get_mass(protein):
    mass_table = {
        "A": 71.03711,
        "C": 103.00919,
        "D": 115.02694,
        "E": 129.04259,
        "F": 147.06841,
        "G": 57.02146,
        "H": 137.05891,
        "I": 113.08406,
        "K": 128.09496,
        "L": 113.08406,
        "M": 131.04049,
        "N": 114.04293,
        "P": 97.05276,
        "Q": 128.05858,
        "R": 156.10111,
        "S": 87.03203,
        "T": 101.04768,
        "V": 99.06841,
        "W": 186.07931,
        "Y": 163.06333
    }
    mass = 0
    for aa in protein:
        mass += mass_table[aa]
    mass = "%.3f" % mass
    return mass


sample = "SKADYEK"
print(get_mass(sample))

with open("rosalind_prtm.txt", "r") as f:
    string = f.read().strip()
    print(get_mass(string))

Output

821.392
115879.301

3、Locating Restriction Sites

Problem

限制性核酸内切酶是可以识别并附着特定的脱氧核苷酸序列，并在每条链中特定部位的两个脱氧核糖核苷酸之间的磷酸二酯键进行切割的一类酶，简称限制酶。限制酶识别DNA序列中的回文序列，回文序列是一种旋转对称结构，在轴的两侧序列相同而反向，特点是在该段的碱基序列的互补链之间正读反读都相同。例如，GCATGC就是一个回文序列，因为它的反向补码也是GCATGC。

Given： FASTA格式的长度最多为1 kbp的DNA字符串。

Return： 在长度为4和12之间的字符串中，每个反向回文的位置和长度。可以按任何顺序输出。

Sample Dataset

>Rosalind_24
TCAATGCATGCGGGTCTATATGCAT

Sample Output

4 6
5 4
6 6
7 4
17 4
18 4
20 6
21 4

example

从第四个碱基开始，长度为6的序列为回文序列。其中注意DNA的互补链序列，取互补碱基之后，需要再反转一下。
**原始：**TCA ATGCAT GCGGGTCTATATGCAT
互补： TACGTA
最终： ATGCAT

注意：回文序列与计算机中的回文字符串并不一样！
回文字符串是一个正读和反读都一样的字符串，比如“level”或者“noon”等等就是回文字符串。

Code

# Locating Restriction Sites

def read_fasta(file):
    """读取数据"""
    sequences = {}
    with open(file, "r") as f:
        for line in f.readlines():
            line = line.strip()
            if line[0] == ">":
                name = line[1:]
                sequences[name] = ""
            else:
                sequences[name] += line
    return sequences


def get_other_seq(sequence):
    """获得DNA互补序列"""
    sequence = sequence[::-1]
    other_seq = ""
    for base in sequence:
        if base == "A":
            other_seq += "T"
        elif base == "T":
            other_seq += "A"
        elif base == "C":
            other_seq += "G"
        elif base == "G":
            other_seq += "C"
    return other_seq


def get_result(seq):
    for i in range(len(seq)):
        # 限制回文字符串的长短
        for j in range(4, 13):
            if i + j <= len(seq):
                cut = seq[i:i + j]
                if cut == get_other_seq(cut):
                    # print(cut)
                    # print(get_other_seq(cut))
                    print(i + 1, j)


seqs = read_fasta("rosalind_revp.txt").values()
seqs = list(seqs)[0]
get_result(seqs)

Output

1 4
1 6
2 4
2 6
3 4
4 4
...
890 6
891 4
898 4
909 4