生信刷题之ROSALIND——Part 4 (MPRT, MRNA, ORF)

最新推荐文章于 2023-08-28 18:41:56 发布

Dzfly..

最新推荐文章于 2023-08-28 18:41:56 发布

阅读量1.2k

点赞数

分类专栏：生信学习文章标签： python ROSALIND 生信

本文链接：https://blog.csdn.net/narutodzx/article/details/130542438

版权

生信学习专栏收录该内容

30 篇文章 144 订阅

订阅专栏

写在前面

本来打算每周更5题，但是题目越来越复杂，想要搞明白，花费的时间越来越多，每个推文的阅读时间也会增长，所以之后会减少题目的数量，同时会在标题标注题目的ID，方便之后复习查找。

1、Finding a Protein Motif

Problem

Motif是一段典型的序列或者一个结构。一般来说，我们称为基序。一般情况下是指构成任何一种特征序列的基本结构。通俗来讲，即是有特征的短序列，一般认为它是拥有生物学功能的保守序列，可能包含特异性的结合位点，或者是涉及某一个特定生物学过程的有共性的序列区段。比如蛋白质的序列特异性结合位点，如核酸酶和转录因子。
为了考虑到蛋白质基序的不同形式，一个蛋白质基序可以用如下的简写来表示:[XY]表示"X或Y"，而{X}表示"除X之外的任何氨基酸"。例如，N-糖基化基序被写成N{P}[ST]{P}。

蛋白质序列可以在UniProt中查找，在https://www.uniprot.org/ 中输入蛋白质access ID（例如：B5ZC00），点击搜索，在sequence中可以看到该蛋白的序列，点击download即可获得fasta格式的蛋白序列信息。具体操作如下图。

fasta信息网址如下，访问其他蛋白质，只需要将名字替换即可。

https://rest.uniprot.org/uniprotkb/{access ID}.fasta

Given： 最多15个UniProt蛋白质数据库access IDs

Return： 对于每个具有N-糖基化基序的蛋白质，输出其给定的access ID，然后输出在该蛋白质中可以找到基序的位置列表。

Sample Dataset

A2Z669
B5ZC00
P07204_TRBM_HUMAN
P20840_SAG1_YEAST

Sample Output

B5ZC00
85 118 142 306 395
P07204_TRBM_HUMAN
47 115 116 382 409
P20840_SAG1_YEAST
79 109 135 248 306 348 364 402 485 501 614

Code

import urllib.request
import re


def get_sequence(accession):
    url = f"https://rest.uniprot.org/uniprotkb/{accession}.fasta"
    request = urllib.request.Request(url)
    response = urllib.request.urlopen(request).read()
    response = str(response)
    # print(response)
    sequence = "".join(response.split("\\n")[1:-1])
    return sequence


def find_motifs(sequence):
    # 正则表达式
    pattern = re.compile(r"N[^P][ST][^P]")
    motif_positions = ""
    # 使用re.finditer有可能会丢失几个位点，感兴趣的可以试一下
    # for match in re.finditer(pattern, sequence):
    #     motif_positions += f"{str(match.start() + 1)} "
    index = 0
    while index < len(sequence):
        index += 1
        # 将序列切割，每次减少一个氨基酸，保证每个位点都被筛选
        if re.search(pattern, sequence[index:]) is None:
            break
        elif re.match(pattern, sequence[index:]) is not None:
            motif_positions += f"{index + 1} "
    return motif_positions


def get_result(file):
    for accession in file.readlines():
        accession = accession.strip()
        # print(accession[:6])
        sequence = get_sequence(accession[:6])
        position = find_motifs(sequence)
        if position:
            print(f"{accession}\n{position}")


with open("code1_example.txt", "r") as f:
    get_result(f)

with open("rosalind_mprt.txt", "r") as f1:
    get_result(f1)

Output

B5ZC00
85 118 142 306 395 
P07204_TRBM_HUMAN
47 115 116 382 409 
P20840_SAG1_YEAST
79 109 135 248 306 348 364 402 485 501 614 
----------
P02974_FMM1_NEIGO
67 68 121 
P01044_KNH1_BOVIN
47 87 168 169 197 204 
Q9QSP4
196 250 326 443 
P02790_HEMO_HUMAN
64 187 240 246 453 
P40225_TPO_HUMAN
197 206 234 255 340 348 
P24592_IBP6_HUMAN
229 
B3ET80
6 
Q55AB5
6 
P01042_KNH_HUMAN
48 169 205 294 
Q1E9Q9
185 255 347 640 1326 
Q9D9T0
154

2、Inferring mRNA from Protein

Problem

翻译（Translation），是蛋白质生物合成（基因表达中的一部分，基因表达还包括转录）过程中的第一步。翻译是根据遗传密码的中心法则，将成熟的mRNA分子（由DNA通过转录而生成）中“碱基的排列顺序”（核苷酸序列）解码，并生成对应的特定氨基酸序列的过程。

当研究人员发现一种新的蛋白质时，他们希望推断出这种蛋白质可能被翻译的mRNA链，从而使他们能够在基因组中定位与这种蛋白质相关的基因。任何RNA链都可以被翻译成一个独特的蛋白质串，但给定蛋白质对应的mRNA却并不唯一，这是因为一个氨基酸往往对应多个密码子。

模算数（Modular arithmetic）是一个整数的算术系统，其中数字超过一定值后（称为模）后会“卷回”到较小的数值。十二小时制就是一种模算数。
模除（modulo）得到的是一个数除以另一个数的余数。

蛋白质过长，对应的mRNA个数会非常多，可以使用模除来避免巨大的数字导致的内存溢出。

Given： 一种长度不超过1000 aa的蛋白质。

Return： 翻译蛋白质的不同RNA链的总数，模值为1,000,000（不要忽视终止密码子在蛋白质翻译中的重要性。）

Sample Dataset

Sample Output

example

例如给定的蛋白质序列为”MA“，M对应的密码子为AUG，A对应的为GCU，GCC，GCA，GCG。且有三个终止密码子(UAA，UAG，UGA)。因此可能的mRNA数量为 1×4×3=12 个。

Code

# Inferring mRNA from Protein

def get_number(sequence):
    rna_codon_table = {
        "F": ["UUU", "UUC"],
        "L": ["UUA", "UUG", "CUU", "CUC", "CUA", "CUG"],
        "I": ["AUU", "AUC", "AUA"],
        "M": ["AUG"],
        "V": ["GUU", "GUC", "GUA", "GUG"],
        "S": ["UCU", "UCC", "UCA", "UCG", "AGU", "AGC"],
        "P": ["CCU", "CCC", "CCA", "CCG"],
        "T": ["ACU", "ACC", "ACA", "ACG"],
        "A": ["GCU", "GCC", "GCA", "GCG"],
        "Y": ["UAU", "UAC"],
        "H": ["CAU", "CAC"],
        "Q": ["CAA", "CAG"],
        "N": ["AAU", "AAC"],
        "K": ["AAA", "AAG"],
        "D": ["GAU", "GAC"],
        "E": ["GAA", "GAG"],
        "C": ["UGU", "UGC"],
        "W": ["UGG"],
        "R": ["CGU", "CGC", "CGA", "CGG", "AGA", "AGG"],
        "G": ["GGU", "GGC", "GGA", "GGG"],
        "Stop": ["UAA", "UAG", "UGA"]
    }
    number = 1
    for aa in sequence:
        if aa in rna_codon_table.keys():
            number *= len(rna_codon_table[aa])
    number *= len(rna_codon_table["Stop"])
    return number % 1000000


print(get_number("MA"))

with open("rosalind_mrna.txt", "r") as f:
    seq = f.read().strip()
    print(get_number(seq))

Output

12
96832

如果不对结果进行模除，结果如下：

2002089841612391663582453283078937275745368421066748680475492424910330611727107825022494435217840287274683628870658010094351191940325350642280590952329513016249187110224840917378827525070051420624021814551693291056364952352113586133745933273636168184838678332737838761697830520116547508864704072632020800420739671680173112061530224686727264605719818801245545285674925560796396475061368452414686546954783174783246059660831096832

3、Open Reading Frames

Problem

开放阅读框【Open reading frame(ORF)】是指在给定的阅读框架中，不包含终止密码子的一串序列。这段序列是生物个体的基因组中，可能作为蛋白质编码序列的部分。基因中的ORF包含并位于开始编码与终止编码之间。由于一段DNA或RNA序列有多种不同读取方式，因此可能同时存在许多不同的开放阅读框架。

一段5’-UCUAAAGGUCCA-3’序列。此序列共有3种读取法：
UCU AAA GGU CCA
CUA AAG GUC
UAA AGG UCC
由于UAA为终止密码子，因此第三种读取法不具编译出蛋白质的潜力，故只有前两者为开放阅读框架。

Given： FASTA格式的长度不超过1 kbp的DNA字符串。

Return： 根据ORF得到的蛋白序列，可以以任意顺序给出。

Sample Dataset

>Rosalind_99
AGCCATGTAGCTAACTCAGGTTACATGGGGATGACCCCGCGACTTGGATTAGAGTCTCTTTTGGAATAAGCCTGAATGATCCGAGTAGCATCTCAG

Sample Output

MLLGSFRLIPKETLIQVAGSSPCNLS
M
MGMTPRLGLESLLE
MTPRLGLESLLE

Code

# Open Reading Frames

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 other_single_chain(sequence):
    """获得互补链序列"""
    # 翻译是5'端到3'端，DNA两条链方向相反
    sequence = sequence[::-1]
    other_sequence = ""
    for base in sequence:
        if base == "A":
            other_sequence += "T"
        elif base == "T":
            other_sequence += "A"
        elif base == "C":
            other_sequence += "G"
        elif base == "G":
            other_sequence += "C"
    return other_sequence


def translation(orf):
    """将所选的orf翻译成蛋白质"""
    dna_codon_table = {
        'ATA': 'I', 'ATC': 'I', 'ATT': 'I', 'ATG': 'M',
        'ACA': 'T', 'ACC': 'T', 'ACG': 'T', 'ACT': 'T',
        'AAC': 'N', 'AAT': 'N', 'AAA': 'K', 'AAG': 'K',
        'AGC': 'S', 'AGT': 'S', 'AGA': 'R', 'AGG': 'R',
        'CTA': 'L', 'CTC': 'L', 'CTG': 'L', 'CTT': 'L',
        'CCA': 'P', 'CCC': 'P', 'CCG': 'P', 'CCT': 'P',
        'CAC': 'H', 'CAT': 'H', 'CAA': 'Q', 'CAG': 'Q',
        'CGA': 'R', 'CGC': 'R', 'CGG': 'R', 'CGT': 'R',
        'GTA': 'V', 'GTC': 'V', 'GTG': 'V', 'GTT': 'V',
        'GCA': 'A', 'GCC': 'A', 'GCG': 'A', 'GCT': 'A',
        'GAC': 'D', 'GAT': 'D', 'GAA': 'E', 'GAG': 'E',
        'GGA': 'G', 'GGC': 'G', 'GGG': 'G', 'GGT': 'G',
        'TCA': 'S', 'TCC': 'S', 'TCG': 'S', 'TCT': 'S',
        'TTC': 'F', 'TTT': 'F', 'TTA': 'L', 'TTG': 'L',
        'TAC': 'Y', 'TAT': 'Y', 'TAA': '_', 'TAG': '_',
        'TGC': 'C', 'TGT': 'C', 'TGA': '_', 'TGG': 'W',
    }
    protein = ""
    for i in range(0, len(orf), 3):
        codon = orf[i:i + 3]
        if codon in dna_codon_table:
            if dna_codon_table[codon] != "_":
                protein += dna_codon_table[codon]
            else:
                break
    return protein


def find_orfs(sequence):
    """根据起始密码子和终止密码子选择可能的orf"""
    orfs = []
    start_codon = 'ATG'
    stop_codons = ['TAA', 'TAG', 'TGA']
    for i in range(len(sequence)):
        if sequence[i:i + 3] == start_codon:
            for j in range(i + 3, len(sequence), 3):
                if sequence[j:j + 3] in stop_codons:
                    orf = sequence[i:j]
                    protein = translation(orf)
                    orfs.append(protein)
                    break
    return orfs


def get_result(file):
    """打印结果"""
    seq = read_fasta(file)
    seq = list(seq.values())[0]
    other_seq = other_single_chain(seq)
    orfs = set()
    orfs.update(find_orfs(seq))
    orfs.update(find_orfs(other_seq))
    for orf in orfs:
        print(orf)


get_result("code3_example.txt")
print("----------")
get_result("rosalind_orf.txt")

Output

M
MTPRLGLESLLE
MGMTPRLGLESLLE
MLLGSFRLIPKETLIQVAGSSPCNLS
----------
MLLHGFSCAYWTFLSRTGGVILGVARYLIFCMWHALRDSFS
MAVNNHTVALVTHHAHTCEVCCG
MLYNNRVR
MSHERHRVVINSHSGPILLRGTWSPSVTKAPRKYTNRPFRSIVDYTYSRWASILP
MELRT
MSSYRSALFKGWSLWRPDPFTT
MACYQTSRAPLYPQHTSHV
MA
MVVMAA
MLGLSVSGDRS
MRSPSCRAPRRIPSYSPRRLKRNHLRCLLRRESCPRYVSRYLSTVNKLAFQRCVQLLNTSRAPSVEEGQENESRNACHMQNIR
MFGPSPLNLGACPG
MGSHVPIGLF
MVKY
MRATCRILDSVRHQE
MKGI
MSARARSKVQRRRAEHGVGSSAPMKGI
MSVMSHERHRVVINSHSGPILLRGTWSPSVTKAPRKYTNRPFRSIVDYTYSRWASILP
MPIDTKTCRPRDFRTVMLGLSVSGDRS
MPASPRVLPSLCQPLLVDRE
MIIPLHYS
MDG
MWHALRDSFS
MSAATCRP
MSSSQHSLGARKGWSFEPRPDVHWLGPIPADRQT
MSSLLAAY
MLRIKLPHPIIVKHTDVSQGTLQGSAATGRTWRRQ
MAA
MPSSRFQDRDARFVGQRGSVLTSGRRV
MVPLQATRTITRYPPRCPARRAPHELIACCLLANKQQAMSSSQHSLGARKGWSFEPRPDVHWLGPIPADRQT
MCLLDFSESHRRSYSWCRTLSNILHVARIARFVLLTFFYARGPGSV
M
MQNIR
MRTHVIAFRMLRIKLPHPIIVKHTDVSQGTLQGSAATGRTWRRQ
MSS