import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, KFold
from sklearn.utils import shuffle
from tensorflow.keras.layers import Dense, Conv2D, BatchNormalization, Activation, add, multiply
from tensorflow.keras.layers import AveragePooling2D, Input, Flatten, GlobalAvgPool2D, Reshape
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.models import Model
from tensorflow.keras.regularizers import l2
from tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping
from scipy.stats import pearsonr
def se_block(input_feature, ration=4):
channel = input_feature.shape[-1]
# squeeze: H*W*C 压缩 1*1*C 大小的特征图,有全局视野
se_feature = GlobalAvgPool2D()(input_feature)
se_feature = Reshape((1,1,channel))(se_feature)
# excitation: 对squeeze后的结果做一个非线性变换,得到不同channel的重要性大小
se_feature = Dense(channel//ration, activation='relu', kernel_initializer='he_normal',
use_bias=True, bias_initializer='zeros')(se_feature)
se_feature = Dense(channel, activation='sigmoid', kernel_initializer='he_normal',
use_bias=True, bias_initializer='zeros', name='weight')(se_feature)
# reweight 特征重标定:用excitation后的结果作为权重,乘到输入特征上
se_feature = multiply([input_feature, se_feature])
return se_feature
# build net
def SimpleNet(input_shape, filter_len):
X_input = Input(input_shape)
# 卷积 filter_num=16 conv_size=filter_len*4 stride=1 l2正则化
X = Conv2D(16, kernel_size=(filter_len, 4), strides=1, padding='valid',
kernel_initializer='glorot_normal', kernel_regularizer=l2(1e-4))(X_input)
# 卷积后的 feature map 标准化
X = BatchNormalization(axis=3, name='bn_conv1')(X)
X = Activation('relu')(X)
y = se_block(X)
# 将 se_block之后的 y 和 X add
X = add([X, y])
X = AveragePooling2D(pool_size=(2,1))(X)
y = Flatten()(X)
# 铺平后加上全连接
outputs = Dense(1, activation='linear', kernel_initializer='glorot_normal')(y)
model = Model(inputs=X_input, outputs=outputs) # save model system
print(model.summary())
return model
# load data
def load_data(data_path, label_path):
x = np.load(data_path)
x = np.expand_dims(x, 3)
y = np.array(pd.read_csv(label_path, header=None).values).reshape([x.shape[0], 1])
X, X_test, Y, Y_test = train_test_split(x, y, test_size=0.2, random_state=123456)
print(X.shape, Y.shape) #(784,100,4,1) (784,1)
return X, X_test, Y, Y_test
if __name__=='__main__':
data_path = 'sgRNA_980_100bp.npy'
label_path = 'label.csv'
X, X_test, Y, Y_test = load_data(data_path, label_path)
# 用十折交叉验证
kfold = KFold(n_splits=10, shuffle=True)
train_cv_result = pd.DataFrame(columns=['filter_len', 'train_pcc']) # to save
# train model
for index, filter_len in enumerate(range(5,26,2)):# filter_size: 5*4、 7*4、 9*4、 11*4...25*4
y_test_all = []
y_pred_all = []
for train_index, test_index in kfold.split(X,Y):
x_train, y_train = X[train_index], Y[train_index] # train in trainSet
x_test, y_test = X[test_index], Y[test_index] # test in trainSet
model = SimpleNet(input_shape=(100,4,1), filter_len=filter_len)
model.compile(loss='mse', optimizer=Adam())
early_stopping = EarlyStopping(monitor='val_loss', patience=10, verbose=0)
# train
model.fit(x_train, y_train, epochs=100, shuffle=True, batch_size=64,
validation_split=0.1, verbose=1, callbacks=[early_stopping])
# predict
y_pred = model.predict(x_test)
y_test_all.extend(y_test)
y_pred_all.extend(y_pred)
# pcc
pcc_train, _ = pearsonr(np.array(y_test_all).flatten(),
np.array(y_pred_all).flatten())
# save best
train_cv_result.loc[index] = [filter_len, pcc_train]
print('When filter length is %d, the pcc in trainSet is %.3f'%(filter_len, pcc_train))
best_length = train_cv_result['filter_len'][train_cv_result['train_pcc'].idxmax(axis=0)]
best_length = int(best_length)
model = SimpleNet(input_shape=(100,4,1), filter_len=best_length)
model.compile(loss='mse', optimizer=Adam())
early_stopping = EarlyStopping(monitor='val_loss', patience=10, verbose=0)
model.fit(X, Y, epochs=100, shuffle=True, batch_size=64,
validation_split=0.1, verbose=1, callbacks=[early_stopping])
Y_pred = model.predict(X_test)
pcc_test, _ = pearsonr(Y_test.ravel(), Y_pred.ravel())
print('The pcc value in independent testSet is %.3f' % pcc_test)
# print weights of 16 filters in SE_block
se_layer_model = Model(inputs=model.input,
outputs=model.get_layer('weight').output)
se_output = se_layer_model.predict(X_test)
print("The weights of 16 filters are:")
print(se_output[0])
Keras实现SEnet实验
最新推荐文章于 2024-05-12 17:35:25 发布