机器学习从蛋白序列预测蛋白分类（一）
一，背景与目标：
随着测序技术的快速发展，GenBank等数据库中存储了大量基因、蛋白序列信息，其中大部分尚无标注，如何充分利用GenBank等数据库现有数据资源，挖掘数据信息，为精准医疗、药物研发等生物大健康领域提供有价值的信息具有重要意义。目前这方面已有BLAST等生物信息技术可用，这里希望尝试机器学习技术在这方面的应用。
二，数据清洗与处理
数据来源于Kaggle竞赛（https://www.kaggle.com/abharg16/predicting-protein-classification）。
1，载入数据
导入需要用到的模块pandas、Numpy及sklearn的相关模块
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score,confusion_matrix,classification_report
pd.set_option('display.max_columns', None) #输出所有列
读入数据
#1) import datasets
df_seq=pd.read_csv('pdb_data_seq.csv')
df_char=pd.read_csv('pdb_data_no_dups.csv')
print('Datasets have been loaded...')
2，查看数据的基本信息
（1）df.head()默认显示数据集的前5行信息
#display data info
print(df_seq.head())
输出如下：
df_seq数据集包含‘StructureId’,‘chainId’,‘sequence’,‘residueCount’,'macromoleculeType’共5列
print(df_char.head())
输出如下：
df_char包含‘structureId’,‘classification’,‘experimentTechnique’,‘macromoleculeType’,‘residueCount’,
‘resolution’,‘structureMoleculeWeight’,‘crystallizatonMethod’,‘crystallizationTempK’,‘densityMattews’,
‘densityPercentSol’，‘pdbxDetails’,‘phValue’,‘publicationYear’’'共14列信息
从‘macromoleculeType’这一列的数据可以看到，数据集包含非蛋白的序列
（2）用df.info()查看数据类型、数据缺失情况
print(df_seq.info())
输出如下：
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 467304 entries, 0 to 467303
Data columns (total 5 columns):
structureId 467304 non-null object
chainId 467294 non-null object
sequence 467276 non-null object
residueCount 467304 non-null int64
macromoleculeType 432487 non-null object
dtypes: int64(1), object(4)
memory usage: 17.8+ MB
None
这里可以看到有几点信息：
首先，df_seq数据集有467304条数据，每条数据有5个维度的信息；
第二，只有‘residueCount’一列为数值型数据（int64）；
第三，只有‘structureId’和‘residueCount’两列无数据缺失(non-null)
print(df_char.info())
输出：
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 141401 entries, 0 to 141400
Data columns (total 14 columns):
structureId 141401 non-null object
classification 141399 non-null object
experimentalTechnique 141401 non-null object
macromoleculeType 137636 non-null object
residueCount 141401 non-null int64
resolution 128589 non-null float64
structureMolecularWeight 141401 non-null float64
crystallizationMethod 96242 non-null object
crystallizationTempK 97039 non-null float64
densityMatthews 124724 non-null float64
densityPercentSol 124749 non-null float64
pdbxDetails 118534 non-null object
phValue 105110 non-null float64
publicationYear 117602 non-null float64
dtypes: float64(7), int64(1), object(6)
memory usage: 15.1+ MB
None
df_char数据集有141401条数据，每条数据有14个维度的信息
3，数据提取与合并
由于df_char数据集中包含非蛋白序列，本项目仅考察蛋白序列的分类，因此需要对数据做抽提与合并
（1）从两个文件（df_seq和df_char）抽提出蛋白序列
#filter for only protein macromoleculeType
protein_chair=df_char[df_char.macromoleculeType=='Protein']
protein_seq=df_seq[df_seq.macromoleculeType=='Protein']
print(protein_seq.head())
输出看一下：
structureId chainId sequence \
4 101M A MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...
7 102L A MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE...
8 102M A MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...
11 103L A MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNSLDAAK...
12 103M A MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...
residueCount macromoleculeType
4 154 Protein
7 165 Protein
8 154 Protein
11 167 Protein
12 154 Protein
（2）从两个表格里筛选需要的列进行数据合并
protein_seq文件和protein_char文件的‘structureId’是唯一的，所以可以作为合并的连接键。
因为我们需要sequence和对应的classification信息，而这两个信息分别处于protein_seq和protein_char两个文件中，所以需要将这两列合并到一个表格，属于表格横向合并，用pd.mergy()函数或者pd.join()函数，先从protein_seq和protein_char中提取出对应的sequence和classification
#select only necessary variables to join
protein_char=protein_char[['structureId','classification']]
protein_seq=protein_seq[['structureId','sequence']]
print(protein_seq.head(2))
print(protein_char.head(2))
输出：
structureId sequence
4 101M MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...
7 102L MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE..
structureId classification
2 101M OXYGEN TRANSPORT
4 102L HYDROLASE(O-GLYCOSYL)
以‘structureId’为连接键用pd.mergy()合并数据
或者先将‘structureId’设为index,再用join()合并
# Join two datasets on structureId using pd.mergy
model_df=protein_seq.set_index('structureId').join(protein_char.set_index('structureId'))
print(model_df.head(2))
输出：
sequence \
structureId
101M MVLSEGEWQLVLHVWAKVEADVAGHGQDILIRLFKSHPETLEKFDR...
102L MNIFEMLRIDEGLRLKIYKDTEGYYTIGIGHLLTKSPSLNAAAKSE...
classification
structureId
101M OXYGEN TRANSPORT
102L HYDROLASE(O-GLYCOSYL)
print(model_df.info())
输出看一下，合并后的表格共有346325条数据，每条数据2个维度，sequence和classification独有数据缺失
<class 'pandas.core.frame.DataFrame'>
Int64Index: 346325 entries, 0 to 346324
Data columns (total 2 columns):
structureId 346325 non-null object
sequence 346322 non-null object
classification 346324 non-null objec
dtypes: object(3)
memory usage: 10.6+ MB
None
（3）查看合并后数据缺失情况，并处理数据缺失:
#check NA
print(model_df.isnull().sum())
输出：
structureId 0
sequence 3
classification 1
dtype: int64
删除缺失数据所在行
model_df=model_df.dropna()
print(model_df.shape[0]) #查看删除后还剩多少条数据
输出：346321，原来有346325条数据，删除4条数据
（4）统计classification分类信息
counts=model_df['classification'].value_counts()
print(counts)
输出
HYDROLASE 46336
TRANSFERASE 36424
OXIDOREDUCTASE 34321
IMMUNE SYSTEM 15615
...
REPLICATION,OXIDOREDUCTASE 1
Virulence Factor 1
Name: classification, Length: 4468, dtype: int64
选择数量超过1000的数据构建模型
print(counts[(counts > 1000)])
输出的是超过1000的classification类型：
HYDROLASE 46336
TRANSFERASE 36424
...
TRANSLATION 1007
BIOSYNTHETIC PROTEIN 1006
Name: classification, dtype: int64
找出classification中统计的数量超过1000的分类名称
print(counts[(counts > 1000)].index)
输出：
Index(['HYDROLASE', 'TRANSFERASE', 'OXIDOREDUCTASE', 'IMMUNE SYSTEM', 'LYASE',
'HYDROLASE/HYDROLASE INHIBITOR', 'TRANSCRIPTION', 'VIRAL PROTEIN',
'TRANSPORT PROTEIN', 'VIRUS', 'SIGNALING PROTEIN', 'ISOMERASE'
数量超过1000的分类名称存储在types中，且数据类型转换为ndarray（即numpy数组类型）
types=np.asarray(counts[(counts > 1000)].index)
查看types中的分类名称是否在‘classification’所在列，isin()返回对应index的True或False
print(model_df.classification.isin(types))
输出：
structureId
101M False
102L False
...
117E True
117E True
118L False
从model_df中抽取满足条件的数据（数量超1000）
data=model_df[model_df.classification.isin(types)]
print(data.head(2))
输出看一下是什么样的
sequence \
structureId
10GS PPYTVVYFPVRGRCAALRMLLADQGQSWKEEVVTVETWQEGSLKAS...
10GS PPYTVVYFPVRGRCAALRMLLADQGQSWKEEVVTVETWQEGSLKAS...
classification
structureId
10GS TRANSFERASE/TRANSFERASE INHIBITOR
10GS TRANSFERASE/TRANSFERASE INHIBITOR
至此，最终用于模型构建的数据集已经抽提结束，看一下一共多少数据
print(data.shape[0])
输出：278866
未完待续…
以下内容见机器学习从蛋白序列预测蛋白分类（二）
三，特征提取
四，模型构建
五，模型评估