Overview1.加载数据2.处理数据3.构建模型4.编译并训练模型 Overview 本篇我们来记录一下怎么使用TensorFlow2.0当中的Keras模块来进行CNN图片分类。 1.加载数据 我们用经典的猫狗分类数据集来做这次图片分类。 import os import numpy as np import matplotlib.pyplot as plt import tensorflow as tf from tensorflow.keras.models import Sequential from tensorflow.keras.layers import Dense, Conv2D, Flatten, Dropout, MaxPooling2D from tensorflow.keras.preprocessing.image import ImageDataGenerator # 下载数据 _URL = 'https://storage.googleapis.com/mledu-datasets/cats_and_dogs_filtered.zip' path_to_zip = tf.keras.utils.get_file('cats_and_dogs.zip', origin=_URL, extract=True) PATH = os.path.join(os.path.dirname(path_to_zip), 'cats_and_dogs_filtered') 统计下载到本地的数据： train_dir = os.path.join(PATH, 'train') validation_dir = os.path.join(PATH, 'validation') train_cats_dir = os.path.join(train_dir, 'cats') # 训练集猫图 train_dogs_dir = os.path.join(train_dir, 'dogs') # 训练集狗图 validation_cats_dir = os.path.join(validation_dir, 'cats') # 验证集猫图 validation_dogs_dir = os.path.join(validation_dir, 'dogs') # 验证集狗图 num_cats_train = len(os.listdir(train_cats_dir)) num_dogs_train = len(os.listdir(train_dogs_dir)) num_cats_val = len(os.listdir(validation_cats_dir)) num_dogs_val = len(os.listdir(validation_dogs_dir)) total_train = num_cats_train + num_dogs_train total_val = num_cats_val + num_dogs_val 2.处理数据 如果直接拿训练集去训练，会遇到过拟合问题，所以我们需要对训练集进行数据增强。图片的数据增强是很简单的，有水平或垂直翻转，旋转，缩放等等方式。 batch_size = 128 epochs = 15 IMG_HEIGHT = 150 IMG_WIDTH = 150 train_image_generator = ImageDataGenerator( rescale=1./255, rotation_range=45, width_shift_range=.15, height_shift_range=.15, horizontal_flip=True, zoom_range=0.5 ) # train_image_generator = ImageDataGenerator(rescale=1./255) validation_image_generator = ImageDataGenerator(rescale=1./255) 读取图片： train_data_gen = train_image_generator.flow_from_directory(batch_size=batch_size, directory=train_dir, shuffle=True, target_size=(IMG_HEIGHT, IMG_WIDTH), class_mode='binary') val_data_gen = validation_image_generator.flow_from_directory(batch_size=batch_size, directory=validation_dir, target_size=(IMG_HEIGHT, IMG_WIDTH), class_mode='binary') 3.构建模型 我们用tensorflow.keras.models.Sequential来顺序堆叠卷积神经网络，同时加上Dropout层来防止过拟合。 model = Sequential([ Conv2D(16, 3, padding='same', activation='relu', input_shape=(IMG_HEIGHT, IMG_WIDTH ,3)), MaxPooling2D(), Conv2D(32, 3, padding='same', activation='relu'), MaxPooling2D(), Conv2D(64, 3, padding='same', activation='relu'), MaxPooling2D(), Dropout(0.2), # dropout层是用来防止过拟合的 Flatten(), Dense(512, activation='relu'), Dense(1) ]) 查看一下网络结构： model.summary() 显示如下： Model: "sequential" _________________________________________________________________ Layer (type) Output Shape Param # ================================================================= conv2d (Conv2D) (None, 150, 150, 16) 448 _________________________________________________________________ max_pooling2d (MaxPooling2D) (None, 75, 75, 16) 0 _________________________________________________________________ conv2d_1 (Conv2D) (None, 75, 75, 32) 4640 _________________________________________________________________ max_pooling2d_1 (MaxPooling2 (None, 37, 37, 32) 0 _________________________________________________________________ conv2d_2 (Conv2D) (None, 37, 37, 64) 18496 _________________________________________________________________ max_pooling2d_2 (MaxPooling2 (None, 18, 18, 64) 0 _________________________________________________________________ dropout (Dropout) (None, 18, 18, 64) 0 _________________________________________________________________ flatten (Flatten) (None, 20736) 0 _________________________________________________________________ dense (Dense) (None, 512) 10617344 _________________________________________________________________ dense_1 (Dense) (None, 1) 513 ================================================================= Total params: 10,641,441 Trainable params: 10,641,441 Non-trainable params: 0 4.编译并训练模型 model.compile(optimizer='adam', loss=tf.keras.losses.BinaryCrossentropy(from_logits=True), metrics=['acc']) history = model.fit( train_data_gen, steps_per_epoch=total_train // batch_size, epochs=epochs, validation_data=val_data_gen, validation_steps=total_val // batch_size ) 这里，选择优化算法时，一般情况下，无脑选择adam即可。 如果运行训练时报如下错误： ImportError: Could not import PIL.Image. The use of `load_img` requires PIL. 那么安装pillow包即可： pip install pillow 如果运行正常，显示如下： Epoch 1/15 15/15 [==============================] - 31s 2s/step - loss: 0.9137 - acc: 0.4995 - val_loss: 0.6909 - val_acc: 0.5011 Epoch 2/15 15/15 [==============================] - 30s 2s/step - loss: 0.6883 - acc: 0.5016 - val_loss: 0.6954 - val_acc: 0.6205 Epoch 3/15 15/15 [==============================] - 30s 2s/step - loss: 0.6933 - acc: 0.5123 - val_loss: 0.6849 - val_acc: 0.4911 Epoch 4/15 15/15 [==============================] - 30s 2s/step - loss: 0.6825 - acc: 0.5037 - val_loss: 0.6637 - val_acc: 0.6105 Epoch 5/15 15/15 [==============================] - 30s 2s/step - loss: 0.6613 - acc: 0.5716 - val_loss: 0.6391 - val_acc: 0.6529 Epoch 6/15 15/15 [==============================] - 30s 2s/step - loss: 0.6439 - acc: 0.5946 - val_loss: 0.6150 - val_acc: 0.6641 Epoch 7/15 15/15 [==============================] - 31s 2s/step - loss: 0.6196 - acc: 0.6167 - val_loss: 0.6006 - val_acc: 0.6641 Epoch 8/15 15/15 [==============================] - 31s 2s/step - loss: 0.6096 - acc: 0.6282 - val_loss: 0.6135 - val_acc: 0.6708 Epoch 9/15 15/15 [==============================] - 31s 2s/step - loss: 0.6063 - acc: 0.6325 - val_loss: 0.5719 - val_acc: 0.6908 Epoch 10/15 15/15 [==============================] - 30s 2s/step - loss: 0.6040 - acc: 0.6485 - val_loss: 0.6070 - val_acc: 0.6920 Epoch 11/15 15/15 [==============================] - 30s 2s/step - loss: 0.5804 - acc: 0.6688 - val_loss: 0.5853 - val_acc: 0.6998 Epoch 12/15 15/15 [==============================] - 30s 2s/step - loss: 0.5908 - acc: 0.6683 - val_loss: 0.5668 - val_acc: 0.6362 Epoch 13/15 15/15 [==============================] - 30s 2s/step - loss: 0.5828 - acc: 0.6667 - val_loss: 0.5698 - val_acc: 0.7054 Epoch 14/15 15/15 [==============================] - 31s 2s/step - loss: 0.5673 - acc: 0.6976 - val_loss: 0.5514 - val_acc: 0.7054 Epoch 15/15 15/15 [==============================] - 30s 2s/step - loss: 0.5627 - acc: 0.6971 - val_loss: 0.5685 - val_acc: 0.7154 查看模型训练过程： acc = history.history['acc'] val_acc = history.history['val_acc'] loss=history.history['loss'] val_loss=history.history['val_loss'] epochs_range = range(epochs) plt.figure(figsize=(8, 8)) plt.subplot(1, 2, 1) plt.plot(epochs_range, acc, label='Training Accuracy') plt.plot(epochs_range, val_acc, label='Validation Accuracy') plt.legend(loc='lower right') plt.title('Training and Validation Accuracy') plt.subplot(1, 2, 2) plt.plot(epochs_range, loss, label='Training Loss') plt.plot(epochs_range, val_loss, label='Validation Loss') plt.legend(loc='upper right') plt.title('Training and Validation Loss') plt.show() 可以看到并没有过拟合，训练还是比较满意的。 本篇文章主要参考了TensorFlow2.0官方文档，Image classification。
Overview问题描述解决方案参考文献 Overview 今天在服务器（Ubuntu 16.04）使用pip安装深度学习包Keras的时候，发现安装结束之后下面有一行提示小字，大约就是pip有新版本可以更新，使用XXX命令更新，就顺带更新了一下，更新完了pip就不能用了... 问题描述 及时是查看下pip的版本 sudo pip --version 也会报下面的错误： Traceback (most recent call last): File "/usr/bin/pip", line 9, in <module> from pip import main ImportError: cannot import name main 解决方案 差了很多资料，发现是pip新版本的一个bug，跟平台无关，Windows下也会出现这个问题（可以查看这里ImportError: cannot import name main when running pip --version command in windows7 32 bit）。 有3个解决办法： 直接修复这个bug，算是临时的解决方案。 用vim编辑/usr/bin/pip文件 sudo vim /usr/bin/pip 可以看到下面的内容 from pip import main if __name__ == '__main__': sys.exit(main()) 把这部分内容注释掉，换成下面的内容： from pip import __main__ if __name__ == '__main__': sys.exit(__main__._main()) 重新运行pip，不在报错了。 对Windows下的pip安装位置不太熟悉，所以不知道应该在哪里修改这个文件。如果不知道在哪里修改，也可以考虑用python命令讲pip降级到没有bug的版本： python -m pip uninstall pip python -m pip install pip==9.0.3 这个没试过，不知道行不行。 也有人说这个bug是由pip的调用方式引起的，在调用pip时，别用： pip install <package> 而是改用（python2）： python -m pip install <package> 或者（python3）： python3 -m pip install <package> 参考文献 ImportError: cannot import name main when running pip --version command in windows7 32 bit
Overview1.归一化代码2.解释说明 Overview 做T4这个项目时，训练集采用了1502条正负样本数据，生成了大量的特征，由于各特征数据变化范围差异很大，于是我们按照重新缩放法做了数据的归一化，使得各特征的取值范围均控制在(0,1)内。然而在预测独立测试集3371条正负样本时，遇到一些疑问：独立测试集怎么办？如果独立测试集归一化，该怎么归一化，以谁为基准归一化，即它的最大值和最小值从何而来？ 首先，毋庸置疑，独立测试集必须经过归一化。假如不归一化，某些特征值将会远大于1，或者远小于，这样样本会被全部判断为正或者全部为负，这点根据SVM和决策树的原理很容易理解。所以，我们查阅了很多资料，发现提到独立测试集归一化的，几乎没有。最终只在stackoverflow上找到一种思路：保存训练集每种特征原始的极值，应用于独立测试集。这样虽然仍会有部分数据“越界”，但是绝大多数数据将会趋于正常。下面将给出python代码，以记录这次工作。编码部分参考了这篇博客。 1.归一化代码 #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 #Copyright Chris & Young def usage(): print "normalize_test.py usage:" print "python normalize_test.py <options> " print "-i,--input: input an independent test dataset feature file." print "-t,--train: input training dataset feature file." print "-o,--ouput: output a normalized independent test dataset feature file." print "-H,--header: input 'T' if the input_file contains header row, else input 'F'." print "-c,--class: input 'T' if the input_file contains class label column, else input 'F'." print "-h,--help: show the help information." import re import pandas as pd import numpy as np import fileinput import sys, getopt from os import listdir from os.path import isfile, join from sklearn import preprocessing opts, args = getopt.getopt(sys.argv[1:], 'i:o:t:H:c:h',['input=','output=','train=','header=','class=','help']) input_file="" output_file="" header_input="" class_input="" train_file="" #col_start这个变量，是归一化开始列的index。根据是否有class列而变化。 col_start=0 for opt, arg in opts: if opt in ('-i','--input'): input_file = arg elif opt in ('-o','--output'): output_file = arg elif opt in ('-t','--train'): train_file=arg elif opt in('-H','--header'): header_input = arg if header_input == "T": header_input='infer' elif header_input == "F": header_input=None elif opt in ('-c','--class'): class_input = arg if class_input=="T": col_start+=1 elif opt in ('-h', '--help'): usage() sys.exit(2) else: usage() sys.exit(2) #读取数据到dataframe中 df_test = pd.read_csv(input_file,header=None) df_train = pd.read_csv(train_file,header=header_input) #将数据框分成两部分：class列单独拿出来，其余部分归一化处理 dataframe_train2 = df_train.ix[:,col_start:] #数据框需要转化为float类型 dataframe_train2 = np.float64(dataframe_train2) dataframe_test = np.float64(df_test) #数据框需要转化为矩阵 matrix_train2=np.array(dataframe_train2) matrix_test=np.array(dataframe_test) #如果两个矩阵列数不同，程序终止退出 [row_train,col_train]=np.shape(matrix_train2) [row_test,col_test]=np.shape(matrix_test) if col_train!=col_test: print "Training or test dataset error!" sys.exit(2) #下面代码是核心归一化代码 min_max_scaler=preprocessing.MinMaxScaler() minmax_train=min_max_scaler.fit_transform(matrix_train2) minmax_test=min_max_scaler.transform(matrix_test) ''' 此注释中代码为测试代码，按列观察训练集和测试集的最大最小值，从而确定是否正确将测试集归一化 train_max=minmax_train.max(axis=0) train_min=minmax_train.min(axis=0) test_max=minmax_test.max(axis=0) test_min=minmax_test.min(axis=0) print train_max,train_min,test_max,test_min ''' #矩阵重新变为数据框 df_final=pd.DataFrame(minmax_test) #将数据框存储为csv格式的文件 df_final.to_csv(output_file,header=None,index=False) 2.解释说明 训练集输入时，程序会按列保存最大值与最小值，只需要sklearn.preprocessing.MinMaxScaler()这个方法，生成一个对象，把它命名为min_max_scaler，这个对象可以提供fit_transform()方法，传入参数为训练集矩阵matrix_train2，这样，这个min_max_scaler对象就可以保存训练集的最大值和最小值，继续使用min_max_scaler的transform()方法，传入matrix_test测试集矩阵，这里要求测试集是裸矩阵，不含header的。这样就可以不用自己保存最大值和最小值矩阵，不用自定义归一化方法了。 这步以后，得到的归一化后的测试集矩阵，即使会有大于1或者小于的元素，但是绝大部分都会在(0,1)内，这样就不会再出现全部的正样本或者负样本的预测结果了。
Overview1. getoptions.py Overview python程序也可以拥有C语言风格的命令行参数，例如-h/--help输出帮助信息，-i/--input提示输入文件，-o/--output提示输出文件等等。除了清晰明了，还有一个好处就是输入参数的顺序可以随意一些。 在python中，实现这个功能的函数跟C/C++中一样，也是getopt().下面用实际的代码展示一下。 1. getoptions.py 首先，我们定义一个usage()函数，用以展示帮助信息： def usage(): print "getoptions.py usage:" print "python getoptions.py <options> " print "-i,--input: the input filename." print "-o,--ouput: the output filename." print "-h,--help: show the help information." 主函数部分如下： #引入sys和getopt包 import sys, getopt opts, args = getopt.getopt(sys.argv[1:],'i:o:h',['input=','output=','help']) inputFile="" outputFile="" for opt, arg in opts: if opt in ('-i','--input'): inputFile = arg elif opt in ('-o','--output'): outputFile = arg elif opt in ('-h', '--help'): usage() sys.exit(2) else: usage() sys.exit(2) ''' 业务逻辑代码 ''' 主函数代码中，我们可以看到 (1). getopt这个函数是在getopt这个包中的。 这个函数接收2或者3个参数：第一个参数是命令行中所有的输入；第二个参数是所有的短option；第三个是可选项，长option。这个可选的标志就是方括号[]. (2). 注意到，i/o后面均有:，而h后面没有:。 这表明，需要输入的option必须加:，而且，在长option后面还要加上赋值号=； 不需要输入的option，不可以加:和=。 (3). 在解析option的时候，短的option前加-，长的option前加--。 更多详细内容，请参考官方文档getopt — C-style parser for command line options。
Overview1 python编码1.1 t34pssm.py1.2 pssm_smoothed_head.py1.3 总结 Overview 上一篇文章python计算smoothed PSSM（一）当中，介绍了以当前氨基酸残基为基点，左右取相同数目的序列，然后叠加计算。Chris介绍，这样的算法有特定的用场：蛋白质后修饰。但是，普通的蛋白质序列提取特征就不太适用了：因为窗口值(smoothed window)只能取奇数，而如果有偶数长度的序列片段包含有特征，这种算法就会漏掉。于是决定写一个新的python脚本，把所有特征全部包含进去。 想法很简单：以当前残基为基点，直接向后连续取w_smth个氨基酸，并叠加计算最后存入新矩阵的当前位置。我做了一个循环，将窗口值不大于w_smth的矩阵全部存入一个相同的新矩阵当中，这样特征就全了。 1 python编码 1.1 t34pssm.py 这部分代码跟前面的代码只有一处不同：不用判断窗口值是否为奇数。这部分内容可参考上一篇内容。 1.2 pssm_smoothed_head.py 这部分代码完成的功能如下： 将每条序列的pssm矩阵的左半部分截取，存入矩阵PSSM-orig。 对矩阵PSSM_orig进行叠加操作，生成矩阵PSSM_smth_head_full。 根据需要截取PSSM_smth_head_full的前n个序列，并存入PSSM_smth_head_final。 将PSSM_smth_head_final合并为一行写入文件，每条序列占一行。 代码如下： #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 """ Retrieve smoothed_head PSSM features """ import sys import numpy as np import math import re import fileinput def pssm_smth_head_n(fi,output_smth_head,w_smth_head,n): # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV' w_smth_head=int(w_smth_head) n=int(n) Amino_vec = "ARNDCQEGHILKMFPSTWYV" PSSM = [] PSSM_orig = [] seq_cn = 0 # 读取pssm文件 for line, strin in enumerate(fileinput.input(fi)): if line > 2: str_vec = strin.split()[1:22] if len(str_vec) == 0: break PSSM.append(map(int, str_vec[1:])) seq_cn += 1 print seq_cn fileinput.close() #original PSSM #将每条序列的`pssm`矩阵的左半部分截取，存入矩阵`PSSM-orig` PSSM_orig=np.array(PSSM) #print PSSM_orig PSSM_smth_head_final=np.array([[0.0]*20]*(n*w_smth_head)) #section for PSSM_smth_head features for k in range(1,w_smth_head+1): PSSM_smth_head = np.array([[0.0]*20]*seq_cn) #print PSSM_smth_head #对矩阵`PSSM_orig`进行叠加操作，生成矩阵`PSSM_smth_head_full`。 PSSM_smth_head_full=pssm_smth_head(PSSM_orig,PSSM_smth_head,k,seq_cn) #print PSSM_smth_head_full #print np.shape(PSSM_smth_head_full) #根据需要截取`PSSM_smth_head_full`的前`n`个序列，并存入`PSSM_smth_head_final`。 for i in range(n): PSSM_smth_head_final[i+n*(k-1)]=PSSM_smth_head_full[i] #print PSSM_smth_head_final PSSM_smth_head_final_shp=np.shape(PSSM_smth_head_final) file_out_smth_head=file(output_smth_head,'a') #将`PSSM_smth_head_final`合并为一行写入文件，每条序列占一行 np.savetxt(file_out_smth_head, [np.reshape(PSSM_smth_head_final, (PSSM_smth_head_final_shp[0] * PSSM_smth_head_final_shp[1], ))], delimiter=",") def pssm_smth_head(PSSM_orig,PSSM_smth_head,w_smth_head,l): for i in range(l): if i <=l-w_smth_head: for j in range(i,i+w_smth_head): PSSM_smth_head[i]+=PSSM_orig[j] else: for j in range(i,l): PSSM_smth_head[i]+=PSSM_orig[j] return PSSM_smth_head 1.3 总结 这个程序可以得到比较全的特征。比如取窗口值为10，那么窗口值为1,2,3,4,5,6,7,8,9,10的特征将会全部被包含在内，并合成一行。 在linux中进入文件所在的目录，然后终端运行如下命令： python t34pssm.py T4undrsmp.txt ./t4 ./t4pssm w_smth n 需要保存的结果文件，自己修改。
Overview1 python编码1.1 t34pssm.py1.2 pssm_smoothed.py1.3 总结 Overview 最近几天，Chris和我看了很多论文，对PSSM有了更深的认识。但是，鉴于PSSM本身包含单个位置的信息更明显，而几乎没有包含蛋白质序列片段信息，我们两人思考如何将蛋白质序列片段信息编码，终于找到了一种PSSM的处理方式，这种方式叫做smoothed window,特此记录一下。 该算法原理，请参考这篇论文：Predicting RNA-binding sites of proteins using support vector machines and evolutionary information，并在此感谢该论文作者！并感谢Chris对我的鼓励和帮助！ 1 python编码 1.1 t34pssm.py 这部分代码是主函数，是之前这篇文章中的代码，可以参考这里蛋白质序列特征提取方法之——PSSM，就不详加解释了。 #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 """ python t34pssm.py total_train_60.fasta ./total_train_60 ./total_train_60_pssm w_smth n param: 1. 总的fasta格式蛋白质序列 2. 分开的fasta格式蛋白质序列的文件夹 3. 分开的fasta格式蛋白质序列对应的pssm的文件夹 4.smooth-window 值，要求是奇数。 5.截取的序列长度，一般为25，30,50。本例为30. """ import fileinput import sys from os import listdir from os.path import isfile, join import re from pssm_smoothed import * smplfasta = sys.argv[1] spfasta = sys.argv[2] check_head = re.compile(r'\>') #read from undersample fasta, store smplist = [] smpcnt = 0 for line, strin in enumerate(fileinput.input(smplfasta)): if check_head.match(strin): smplist.append(strin.strip()) smpcnt += 1 onlyfiles = [ f for f in listdir(spfasta) if isfile(join(spfasta,f)) ] fastaDict = {} for fi in onlyfiles: cntnt = '' for line, strin in enumerate(fileinput.input(spfasta+'/'+fi)): if line == 0: cntnt += strin.strip() if cntnt in fastaDict: print strin fastaDict[cntnt] = fi pssmdir = sys.argv[3] w_smth = sys.argv[4] #如果窗口值不是奇数，退出程序 if int(w_smth)%2 ==0: print 'Please change your input argument ' + w_smth + ' to an odd smoothing-window number!!!' sys.exit(1) n=sys.argv[5] for smp in smplist: finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm') for fi in finalist: pssm_single(fi,'total_train_60_pssm_smth',w_smth,n) 1.2 pssm_smoothed.py #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 """ Retrieve smoothed PSSM features """ import sys import numpy as np import math import re import fileinput def pssm_single(fi,output_smth,w_smth,n): # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV' w_smth=int(w_smth) n=int(n) Amino_vec = "ARNDCQEGHILKMFPSTWYV" PSSM = [] PSSM_orig = [] seq_cn = 0 # 读取pssm文件 for line, strin in enumerate(fileinput.input(fi)): if line > 2: str_vec = strin.split()[1:22] if len(str_vec) == 0: break PSSM.append(map(int, str_vec[1:])) seq_cn += 1 print seq_cn fileinput.close() PSSM_smth = np.array([[0.0]*20]*seq_cn) #print PSSM_smth PSSM_orig=np.array(PSSM) #print PSSM_orig #section for PSSM_smth features PSSM_smth_full=pssm_smth(PSSM_orig,PSSM_smth,w_smth,seq_cn) PSSM_smth_final=[[0.0]*20]*n #截取PSSM_smth_full矩阵的前n行，作为输出内容 for i in range(n): PSSM_smth_final[i]=PSSM_smth_full[i] #print PSSM_smth_final[i] PSSM_smth_final_shp=np.shape(PSSM_smth_final) # for i in range(seq_cn): # print PSSM_smth_final[i] file_out_smth=file(output_smth,'a') np.savetxt(file_out_smth, [np.reshape(PSSM_smth_final, (PSSM_smth_final_shp[0] * PSSM_smth_final_shp[1], ))], delimiter=",") #这个函数会求出整条序列的smoothed pssm矩阵 def pssm_smth(PSSM_orig,PSSM_smth,w_smth,l): for i in range(l): #smooth窗口超过pssm上边界 if i <(w_smth-1)/2: for j in range(i+(w_smth-1)/2+1): #print i,j PSSM_smth[i]+=PSSM_orig[j] #print PSSM_smth[i] #smooth窗口超过pssm下边界 elif i>=(l-(w_smth-1)/2): for j in range(i-(w_smth-1)/2,l): #print i,j PSSM_smth[i]+=PSSM_orig[j] #print PSSM_smth[i] else: for j in range(i-(w_smth-1)/2,i+(w_smth-1)/2+1): #print i,j PSSM_smth[i]+=PSSM_orig[j] #print PSSM_smth[i] return PSSM_smth 1.3 总结 该算法以行为单位进行运算。以本行为中心，上下扩展，扩展的上下长度为smooth窗口值，将这些行的值加起来，存入新的矩阵的相同行位置。这样新的矩阵就包含了多个连续氨基酸序列片段的信息，会为特征提取提供新的思路。
Overview1. separatePosNeg.py Overview 机器学习算法在项目中应用的时候，有时候会用到分离fasta格式的正负样本。于是就写了一个python脚本，效果不错，记录一下。 1. separatePosNeg.py #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 import fileinput import sys import re ls=fileinput.input(sys.argv[1]) check_head = re.compile(r'\>') cls = '' smplist = '' for line, strin in enumerate(ls): if check_head.match(strin): strin_vec=strin.split('|') cls = strin_vec[len(strin_vec)-1].strip() if cls == str(sys.argv[2]): smplist+=strin.strip() smplist+='\n' smplist+=ls.next() print smplist 记录下在linux命令行中的用法： #提取正样本时： python separatePosNeg.py original.fasta '1' >> pos.fasta #提取负样本时： python separatePosNeg.py original.fasta '-1' >> neg.fasta
Overview1. t34pssm.py2. pssm_ft.py Overview 昨天跟Chris讨论SVM分类预测准确性的时候，知道PSSM_AC的作用比PSSM作用更明显，于是决定将以前的python脚本改进一下，输出PSSM和PSSM_AC这两个文件，方便观察。该脚本包括两部分，本文将按顺序记录下来。 以前的脚本可以参考我之前的文章蛋白质序列特征提取方法之——PSSM。 1. t34pssm.py #! /usr/bin/env python import fileinput import sys from os import listdir from os.path import isfile, join import re from pssm_ft import * smplfasta = sys.argv[1] spfasta = sys.argv[2] check_head = re.compile(r'\>') #read from undersample fasta, store smplist = [] smpcnt = 0 for line, strin in enumerate(fileinput.input(smplfasta)): if check_head.match(strin): smplist.append(strin.strip()) smpcnt += 1 onlyfiles = [ f for f in listdir(spfasta) if isfile(join(spfasta,f)) ] fastaDict = {} for fi in onlyfiles: cntnt = '' for line, strin in enumerate(fileinput.input(spfasta+'/'+fi)): if line == 0: cntnt += strin.strip() if cntnt in fastaDict: print strin fastaDict[cntnt] = fi pssmdir = sys.argv[3] finalist = [] for smp in smplist: finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm') for fi in finalist: pssm(fi, 'total_train_60.pssm','total_train_60.pssm_AC') #此处多了一个参数，多了一个输出文件 2. pssm_ft.py import sys import numpy as np import math import re import fileinput GROUP = 10 def scale(x, max, min): """ scale x to [-1,1] """ try: return -(1-(x-min)/(max-min))+((x-min)/(max-min)) except Exception: print "ZeroDivisionError, max-min:\s", max-min def pssm_ac_cal(PSSM, g, l): """ :param PSSM: ['1','-1'...] :param g: group :param l: sequence length :return PSSM_AC: g * 20 matrix """ PSSM_AC = np.array([ [0.0] * 20 ] * g) for pg in range(g): l_g = l - pg - 1 for pj in range(20): sum_jl = 0.0 for i in range(l): sum_jl += PSSM[i][pj] sum_jl /= l pssm_acjg = 0.0 for i in range(l_g): pssm_acjg += (PSSM[i][pj]-sum_jl) * (PSSM[i+pg+1][pj]-sum_jl) pssm_acjg /= l_g PSSM_AC[pg][pj] = pssm_acjg return PSSM_AC def pssm(fi,output,output_AC): # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV' Amino_vec = "ARNDCQEGHILKMFPSTWYV" PSSM = [] ACC = [ [0.0] * 20 ] * 20 seq_cn = 0 for line, strin in enumerate(fileinput.input(fi)): if line > 2: str_vec = strin.split()[1:22] if len(str_vec) == 0: break PSSM.append(map(int, str_vec[1:])) seq_cn += 1 #test1_1 += float(str_vec[1]) test1_1 = -103 ACC[Amino_vec.index(str_vec[0])] = map(sum, zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])])) fileinput.close() # section for ACC features ACC_np = np.array(ACC) ACC_np = np.divide(ACC_np, seq_cn) amcnt_max = np.amax(ACC_np) amcnt_min = np.amin(ACC_np) vfunc = np.vectorize(scale) ACC_np = vfunc(ACC_np, amcnt_max,amcnt_min) ACC_shp = np.shape(ACC_np) #section for PSSM_AC features PSSM = np.array(PSSM) PSSM_AC = pssm_ac_cal(PSSM, GROUP, seq_cn) PSSM_shp = np.shape(PSSM_AC) file_out = file(output,'a') file_out_AC = file(output_AC,'a') #此处将保存文件的代码拆分开 np.savetxt(file_out, [np.reshape(ACC_np, (ACC_shp[0] * ACC_shp[1], ))], delimiter=",") np.savetxt(file_out_AC, [np.reshape(PSSM_AC, (PSSM_shp[0] * PSSM_shp[1], ))], delimiter=",")
Overview解决方案 Overview 在转换fasta格式的文件为chen's format文件时，发现前人的程序有些bug，会将最后一条正样本的class标记为-1于是将程序改了一下，这下就没有问题了。 解决方案 fasta格式如下，第一行为头信息，第二行为氨基酸序列： >sp|Q2YIT7|VIRB3_BRUA2 Type IV secretion system protein virB3 OS=Brucella abortus (strain 2308) GN=virB3 PE=3 SV=1|1 MTTAPQESNARSAGYRGDPIFKGCTRPAMLFGVPVIPLVIVGGSIVLLSVWISMFILPLIVPIVLVMRQITQTDDQMFRLLGLKAQFRLIHFNRTGRFWRASAYSPIAFTKRKRES 现在需要将其转换为如下格式，以供cksaap以及aaindex算法之用： MTTAPQESNARSAGYRGDPIFKGCTRPAMLFGVPVIPLVIVGGSIVLLSVWISMFILPLIVPIVLVMRQITQTDDQMFRLLGLKAQFRLIHFNRTGRFWRASAYSPIAFTKRKRES Q2YIT7 any 1 修改后的python代码如下： import sys import re import fileinput check_head = re.compile(r'\>') seqs = '' head_line = '' cls = '' for line, strin in enumerate(fileinput.input(sys.argv[1])): if check_head.match(strin): strin_vec=strin.split('|') cls = strin_vec[len(strin_vec)-1].strip() head_line = strin_vec[1].strip() else: seqs = strin.strip() print seqs+'\t'+head_line+'\t'+'any'+'\t'+cls 有些程序处理类似的格式时，有可能会把下一行或者上一行的信息和本行信息混淆，所以处理完之后最好手动检查一下正负样本交界处的数据。