MIALAB/ACPred-BMF

1# -*- coding: utf-8 -*-2"""3Spyder Editor4 5This is a temporary script file.6"""7 8import pandas as pd9import numpy as np10import os11import warnings12warnings.filterwarnings('ignore')13#######Standardization14from sklearn.preprocessing import StandardScaler15from sklearn import preprocessing16from tensorflow.keras.layers import LSTM,Dense,Activation,Dropout,Bidirectional17from tensorflow.keras import Sequential18from keras.models import *19 20 21def read_fasta(input): #Define a function read_fasta() using def and pass the parameter using the variable input.22    with open(input,'r') as f: # open file23        fasta = [] # Define an empty dictionary24        for line in f:25            line = line.strip() # Remove the trailing newline character.26            if line[0] == '>':27                header = line[1:]28            else:29                sequence = line30                fasta.append(sequence)31    return fasta32##################################################### Read the file using the written function33 34import sys35input_fasta = sys.argv[1] # file type: fasta36output_dir = sys.argv[2]  # file dir37filename = sys.argv[3] #file name38 39 40fa= read_fasta(input_fasta)41 42####Convert to DataFrame43df=pd.DataFrame(pd.Series(fa),columns=['seq'])44df=df.reset_index().rename(columns={'index':'id'})45#df['leixing']=np.where(df['id'].apply(lambda x: '|1'in x),1,0)46df.head()47 48 49feature_csv = os.path.join(os.path.dirname(os.path.abspath(__file__)),'amino_acids_1.csv')50pp_aa=pd.read_csv(feature_csv,encoding='gbk')51pp_aa.head()52embed=pd.get_dummies(pp_aa.iloc[:,1:12])53stdScale=StandardScaler().fit(embed.iloc[:,0:6])54#minmaxScale=preprocessing.MinMaxScaler().fit(embde.iloc[:,0:6])55embed.head()56embed.iloc[:,0:6]=stdScale.transform(embed.iloc[:,0:6])57embed.head()58 59##mixed nature60def encode(x,encode_len,mask):61    seq_encode=list()62    acid=['G','A','V','L','I','F','W','Y','D','H','N','E','K','Q','M','R','S','T','C','P']63    if len(x)>=encode_len:64        for i in range(encode_len):65            m=[0]*36          66            for j in range(20):67                if x[i]==acid[j]:68                    m=embed.iloc[j,:]69            m=list(m)70            seq_encode.append(m)71    else:72        for i in range(len(x)):73            m=[0]*3674            for j in range(20):75                if x[i]==acid[j]:76                    m=embed.iloc[j,:]77            m=list(m)78            seq_encode.append(m)79        for j in range(encode_len-len(x)):80            seq_encode.append([mask]*36)81            #bit_opf.append([0]*30)82    return(seq_encode) 83    84single_n=3685mask_value=10086bb=df['seq']87bb=pd.DataFrame(bb)88#max_num=np.max(bb.seq.apply(lambda x:len(x)))89max_num=5090bb['encode']=bb.seq.apply(lambda x:encode(x,encode_len=max_num,mask=mask_value))91#print(bb)92x_code=list()93for i in range(bb.shape[0]):94    x_code+=bb.iloc[i,1]95 96###Construct data format suitable for LSTM97###Where bb[i, j, :] represents the (j+1)-th amino acid of the (i+1)-th sequence98x=np.array(x_code).reshape(bb.shape[0],max_num,single_n)99#x_test1=np.array(x_code1).reshape(cc.shape[0],max_num,single_n)100#y=df['leixing']101#y_test1=main_ts['leixing']102#x_quanc=x[:,:,0:26]103 104#####################model prediction105##load model106alter_h5 = os.path.join(os.path.dirname(os.path.abspath(__file__)),'alt.h5')107alter_model=load_model(alter_h5)108 109prediction=pd.DataFrame(0,columns = ['index','sequence','prediction'],110                                   index=range(df.shape[0]))111y_pred=alter_model.predict(x)112y_pred_xgb= np.argmax(y_pred, axis=1)113prediction['sequence']=df['seq']114prediction['prediction']=y_pred_xgb115prediction['index'] = prediction.index116 117output_path = output_dir + "/" + filename + "_alter_prediction.csv"118prediction.to_csv(output_path,index=False)119