python 生物信息学数据管理

2021-03-05 本文已影响0人所以suoyi

2021/03/05

一、计算胰岛素序列中的氨基酸频率

任务：计算在这条蛋白质序列上20种氨基酸出现的频率是多少？
思考：1、有多少半胱氨酸可以组成二硫键？
2、是否有特殊的非极性残基数目表明存在可跨膜的结构域？
3、是否有大量阳性电荷残基可以参与核酸结合？

insulin = 'GlVEQCCTSICSLYQLENYCNFVNQHLCGSHLVEALYLVCGERGFFYTPKT'
amino_acid = 'ACDEFGHIKLMNPQRSTVWY'
for aa in amino_acid:
    number = insulin.count(aa)
    print(aa, number)

coden_table = {
    'GCU': 'Ala', 'GCC': 'Ala', 'GCA': 'Ala', 'GCG': 'Ala',
    'CGU': 'Arg', 'CGC': 'Arg', 'CGA': 'Arg', 'CGG': 'Arg', 'AGA': 'Arg', 'AGG': 'Arg',
    'UCU': 'Ser', 'UCC': 'Ser', 'UCA': 'Ser', 'UCG': 'Ser', 'AGU': 'Ser', 'AGC': 'Ser',
    'AUU': 'Ile', 'AUC': 'Ile', 'AUA': 'Ile',
    'UUA': 'Leu', 'UUG': 'Leu', 'CUU': 'Leu', 'CUC': 'Leu', 'CUA': 'Leu', 'CUG': 'Leu',
    'GGU': 'Gly', 'GGC': 'Gly', 'GGA': 'Gly', 'GGG': 'Gly',
    'GUU': 'Val', 'GUC': 'Val', 'GUA': 'Val', 'GUG': 'Val',
    'ACU': 'Thr', 'ACC': 'Thr', 'ACA': 'Thr', 'ACG': 'Thr',
    'CCU': 'Pro', 'CCC': 'Pro', 'CCA': 'Pro', 'CCG': 'Pro',
    'AAU': 'Asn', 'AAC': 'Asn',
    'GAU': 'Asp', 'GAC': 'Asp',
    'UGU': 'Cys', 'UGC': 'Cys',
    'CAA': 'Gln', 'CAG': 'Gln',
    'GAA': 'Glu', 'GAG': 'Glu',
    'CAU': 'His', 'CAC': 'His',
    'AAA': 'Lys', 'AAG': 'Lys',
    'UUU': 'Phe', 'UUC': 'Phe',
    'UAU': 'Tyr', 'UAC': 'Tyr',
    'AUG': 'Met',
    'UGG': 'Trp',
    'UAG': 'Ter', 'UGA': 'Ter', 'UAA': 'Ter'
    }
coden_trans = {
        'Ala': 'A', 'Arg': 'R', 'Asn': 'N', 'Asp': 'D', 'Cys': 'C', 'Gln': 'Q',
        'Glu': 'E', 'Gly': 'G', 'His': 'H', 'Ile': 'I', 'Leu': 'L', 'Lys': 'K', 
        'Met': 'M', 'Phe': 'F', 'Pro': 'P', 'Ser': 'S', 'Thr': 'T', 'Trp': 'W', 
        'Tyr': 'Y', 'Val': 'V', 'Ter': 'X'}

二硫键（disulfide bond）是连接不同肽链或同一肽链中，两个不同半胱氨酸残基的巯基的化学键。二硫键是比较稳定的共价键，在蛋白质分子中，起着稳定肽链空间结构的作用。二硫键数目越多，蛋白质分子对抗外界因素影响的稳定性就愈大。
蛋白质序列中跨越细胞膜的区域，通常为α-螺旋结构，约20~25个氨基酸残基构成他的蛋白质的氨基酸大部分是疏水性氨基酸。
非极性氨基酸（疏水氨基酸）共8种：
丙氨酸（Ala）；缬氨酸（Val）；亮氨酸（Leu）；异亮氨酸（Ile）；
脯氨酸（Pro）；苯丙氨酸（Phe）；色氨酸（Trp）；蛋氨酸（Met）
极性氨基酸（亲水氨基酸）共14种：
极性不带电荷（中性氨基酸）：
甘氨酸（Gly）；丝氨酸（Ser）；苏氨酸（Thr）；半胱氨酸（Cys）；
酪氨酸（Tyr）；天冬酰胺（Asn）；谷氨酰胺（Gln）；
硒半胱氨酸（Sec）；吡咯赖氨酸（Pyl）
极性带正电荷的氨基酸（碱性氨基酸）：
赖氨酸（Lys）；精氨酸（Arg）；组氨酸（His）
极性带负电荷的氨基酸（酸性氨基酸）：
天冬氨酸（Asp）；谷氨酸（Glu) 。

二、整合质谱数据，整合到代谢通路中

蛋白质谱技术简单来说就是一种将质谱仪用于研究蛋白质的技术，通过和理论上蛋白质经过胰蛋白酶消化后产生的一级质谱峰图和二级质谱峰图进行比对而鉴定蛋白质。
质谱实验的最终结果应包含研究样本中检测到的肽(即表达结果)的列表，通过特定数据分析软件(如 Mascot http://www.matrixscience.com/ ) ，质谱测定的肽可以与UniProt 序列匹配，使得输出列表可以用 UniProt ID 的形式表示，并通常存储于CSV（用逗号分隔的)文本文件中。
代谢通路（如Reactome https://reactome.org/）

list1 = ['uniprot id1', 'uniprot id2', 'uniprot id3', 'uniprot id4']
list2 = ['某代谢通路包含的蛋白']
# 此问题是求得两个数据集共同包含的蛋白质
common = []
no_in = []
for i in list1:
    if i in list2:
        common.append(i)
    else:
        no_in.append(i)

三、在蛋白质序列中搜索磷酸化模体

蛋白质磷酸化：指由蛋白质激酶催化的把ATP的磷酸基转移到底物蛋白质氨基酸残基（丝氨酸、苏氨酸、酪氨酸）上的过程，或者在信号作用下结合GTP，是生物体内一种普通的调节方式，在细胞信号转导的过程中起重要作用。蛋白质磷酸化是调节和控制蛋白质活力和功能的最基本、最普遍，也是最重要的机制。

import re
seq = 'VSVLTMFRYAGWLDRLYMLVGTQLAAIIHGVALPLMMLI'
# [ST]Q : 丝氨酸/苏氨酸磷酸化模体
match = re.search('[ST]Q', seq)  # 返回第一个
print(match)  # <re.Match object; span=(21, 23), match='TQ'> 
print(match.start(), match.end())  # 21 23
match = re.findall('[ST]Q', seq)  # ['TQ']
print(match)

seq = 'VSVLTMFRYAGWLDRLYMLVGTQLAAIIHGSQLPLMMLI'
for i in re.finditer('[ST]Q', seq):
    print(i)
# <re.Match object; span=(21, 23), match='TQ'>
# <re.Match object; span=(30, 32), match='SQ'>

四、核糖体的核苷酸频率

创建一个柱状图

23S 核糖体亚基的核酸数目.png

from pylab import figure, title, xlabel, ylabel, xticks, bar, legend, axis, savefig
nucleotides = ['A', 'G', 'C', 'U']
counts = [[606, 1024, 759, 398], [762, 912, 639, 591]]
figure()
title('RNA nucleotides in the ribosome')
xlabel('RNA')  # x轴 变量名
ylabel('base count')  # y轴 变量名
x1 = [2.0, 4.0, 6.0, 8.0]
x2 = [x - 0.5 for x in x1]   # 坐标轴间距
xticks(x1, nucleotides)
bar(x1, counts[1], width=0.5, color="#cccccc", label="E.coli 23S")
bar(x2, counts[0], width=0.5, color="#808080", label="T.thermophilus 23S")
legend()
axis([1.0, 9.0, 0, 1200])
savefig('barplot.png')

barplot.png

五、序列比对

wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.11.0/ncbi-blast-2.11.0+-x64-linux.tar.gz
tar -zxvf ncbi-blast-2.11.0+-x64-linux.tar.gz
echo "export PATH=/home/user/blast/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc
blastn -version

import os
commond = 'blastn  -db hg19.fa  -query example_fasta.txt -evalue 1e-5  
          -perc_identity 100  -task megablast  -out result.blast  
          -outfmt "7 qacc sacc evalue sstart send" '
os.system(commond)

blast_result.png