text: #!/usr/bin/env python
  # coding: utf-8
  
  # # ORF recognition by CNN
  # 
  # Use variable number of bases between START and STOP. Thus, ncRNA will have its STOP out-of-frame or too close to the START, and pcRNA will have its STOP in-frame and far from the START.
  
  # In[1]:
  
  
  import time 
  t = time.time()
  time.strftime('%Y-%m-%d %H:%M:%S %Z', time.localtime(t))
  
  
  # In[2]:
  
  
  PC_SEQUENCES=10000   # how many protein-coding sequences
  NC_SEQUENCES=10000   # how many non-coding sequences
  PC_TESTS=1000
  NC_TESTS=1000
  RNA_LEN=32            # how long is each sequence
  CDS_LEN=16            # min CDS len to be coding
  ALPHABET=4          # how many different letters are possible
  INPUT_SHAPE_2D = (RNA_LEN,ALPHABET,1) # Conv2D needs 3D inputs
  INPUT_SHAPE = (RNA_LEN,ALPHABET) # Conv1D needs 2D inputs
  FILTERS = 16   # how many different patterns the model looks for
  NEURONS = 16
  DROP_RATE = 0.0
  WIDTH = 3   # how wide each pattern is, in bases
  STRIDE_2D = (1,1)  # For Conv2D how far in each direction
  STRIDE = 1 # For Conv1D, how far between pattern matches, in bases
  EPOCHS=100  # how many times to train on all the data
  SPLITS=3  # SPLITS=3 means train on 2/3 and validate on 1/3 
  FOLDS=3  # train the model this many times (range 1 to SPLITS)
  
  
  # In[3]:
  
  
  import sys
  IN_COLAB = False
  try:
      from google.colab import drive
      IN_COLAB = True
  except:
      pass
  if IN_COLAB:
      print("On Google CoLab, mount cloud-local file, get our code from GitHub.")
      PATH='/content/drive/'
      #drive.mount(PATH,force_remount=True)  # hardly ever need this
      #drive.mount(PATH)    # Google will require login credentials
      DATAPATH=PATH+'My Drive/data/'  # must end in "/"
      import requests
      r = requests.get('https://raw.githubusercontent.com/ShepherdCode/Soars2021/master/SimTools/RNA_describe.py')
      with open('RNA_describe.py', 'w') as f:
          f.write(r.text)  
      from RNA_describe import ORF_counter
      from RNA_describe import Random_Base_Oracle
      r = requests.get('https://raw.githubusercontent.com/ShepherdCode/Soars2021/master/SimTools/RNA_prep.py')
      with open('RNA_prep.py', 'w') as f:
          f.write(r.text)  
      from RNA_prep import prepare_inputs_len_x_alphabet
  else:
          print("CoLab not working. On my PC, use relative paths.")
          DATAPATH='data/'  # must end in "/"
          sys.path.append("..") # append parent dir in order to use sibling dirs
          from SimTools.RNA_describe import ORF_counter,Random_Base_Oracle
          from SimTools.RNA_prep import prepare_inputs_len_x_alphabet
  
  MODELPATH="BestModel"  # saved on cloud instance and lost after logout
  #MODELPATH=DATAPATH+MODELPATH  # saved on Google Drive but requires login
  
  
  # In[4]:
  
  
  from os import listdir
  import csv
  from zipfile import ZipFile
  
  import numpy as np
  import pandas as pd
  from scipy import stats  # mode
  
  from sklearn.preprocessing import StandardScaler
  from sklearn.model_selection import KFold
  from sklearn.model_selection import cross_val_score
  
  from keras.models import Sequential
  from keras.layers import Dense,Embedding,Dropout
  from keras.layers import Conv1D,Conv2D
  from keras.layers import Flatten,MaxPooling1D,MaxPooling2D
  from keras.losses import BinaryCrossentropy
  # tf.keras.losses.BinaryCrossentropy
  
  import matplotlib.pyplot as plt
  from matplotlib import colors
  mycmap = colors.ListedColormap(['red','blue'])  # list color for label 0 then 1
  np.set_printoptions(precision=2)
  
  
  # In[5]:
  
  
  rbo=Random_Base_Oracle(RNA_LEN,True)
  pc_all,nc_all = rbo.get_partitioned_sequences(CDS_LEN,10) # just testing
  pc_all,nc_all = rbo.get_partitioned_sequences(CDS_LEN,PC_SEQUENCES+PC_TESTS)
  print("Use",len(pc_all),"PC seqs")
  print("Use",len(nc_all),"NC seqs")
  
  
  # In[6]:
  
  
  # Describe the sequences
  def describe_sequences(list_of_seq):
      oc = ORF_counter()
      num_seq = len(list_of_seq)
      rna_lens = np.zeros(num_seq)
      orf_lens = np.zeros(num_seq)
      for i in range(0,num_seq):
          rna_len = len(list_of_seq[i])
          rna_lens[i] = rna_len
          oc.set_sequence(list_of_seq[i])
          orf_len = oc.get_max_orf_len()
          orf_lens[i] = orf_len
      print ("Average RNA length:",rna_lens.mean())
      print ("Average ORF length:",orf_lens.mean())
      
  print("Simulated sequences prior to adjustment:")
  print("PC seqs")
  describe_sequences(pc_all)
  print("NC seqs")
  describe_sequences(nc_all)
  
  
  # In[7]:
  
  
  pc_train=pc_all[:PC_SEQUENCES]
  nc_train=nc_all[:NC_SEQUENCES]
  pc_test=pc_all[PC_SEQUENCES:]
  nc_test=nc_all[NC_SEQUENCES:]
  
  
  # In[8]:
  
  
  # Use code from our SimTools library.
  X,y = prepare_inputs_len_x_alphabet(pc_train,nc_train,ALPHABET) # shuffles
  print("Data ready.")
  
  
  # In[9]:
  
  
  def make_DNN():
      print("make_DNN")
      print("input shape:",INPUT_SHAPE)
      dnn = Sequential()
      #dnn.add(Embedding(input_dim=INPUT_SHAPE,output_dim=INPUT_SHAPE)) 
      dnn.add(Conv1D(filters=FILTERS,kernel_size=WIDTH,strides=STRIDE,padding="same",
              input_shape=INPUT_SHAPE))
      dnn.add(Conv1D(filters=FILTERS,kernel_size=WIDTH,strides=STRIDE,padding="same"))
      dnn.add(MaxPooling1D())
      dnn.add(Conv1D(filters=FILTERS,kernel_size=WIDTH,strides=STRIDE,padding="same"))
      dnn.add(Conv1D(filters=FILTERS,kernel_size=WIDTH,strides=STRIDE,padding="same"))
      dnn.add(MaxPooling1D())
      dnn.add(Flatten())
      dnn.add(Dense(NEURONS,activation="sigmoid",dtype=np.float32))   
      dnn.add(Dropout(DROP_RATE))
      dnn.add(Dense(1,activation="sigmoid",dtype=np.float32))   
      dnn.compile(optimizer='adam',
                  loss=BinaryCrossentropy(from_logits=False),
                  metrics=['accuracy'])   # add to default metrics=loss
      dnn.build(input_shape=INPUT_SHAPE)
      #ln_rate = tf.keras.optimizers.Adam(learning_rate = LN_RATE)
      #bc=tf.keras.losses.BinaryCrossentropy(from_logits=False)
      #model.compile(loss=bc, optimizer=ln_rate, metrics=["accuracy"])
      return dnn
  model = make_DNN()
  print(model.summary())
  
  
  # In[10]:
  
  
  from keras.callbacks import ModelCheckpoint
  def do_cross_validation(X,y):
      cv_scores = []
      fold=0
      mycallbacks = [ModelCheckpoint(
          filepath=MODELPATH, save_best_only=True, 
          monitor='val_accuracy', mode='max')]   
      splitter = KFold(n_splits=SPLITS)  # this does not shuffle
      for train_index,valid_index in splitter.split(X):
          if fold < FOLDS:
              fold += 1
              X_train=X[train_index] # inputs for training
              y_train=y[train_index] # labels for training
              X_valid=X[valid_index] # inputs for validation
              y_valid=y[valid_index] # labels for validation
              print("MODEL")
              # Call constructor on each CV. Else, continually improves the same model.
              model = model = make_DNN()
              print("FIT")  # model.fit() implements learning
              start_time=time.time()
              history=model.fit(X_train, y_train, 
                      epochs=EPOCHS, 
                      verbose=1,  # ascii art while learning
                      callbacks=mycallbacks,   # called at end of each epoch
                      validation_data=(X_valid,y_valid))
              end_time=time.time()
              elapsed_time=(end_time-start_time)                        
              print("Fold %d, %d epochs, %d sec"%(fold,EPOCHS,elapsed_time))
              # print(history.history.keys())  # all these keys will be shown in figure
              pd.DataFrame(history.history).plot(figsize=(8,5))
              plt.grid(True)
              plt.gca().set_ylim(0,1) # any losses > 1 will be off the scale
              plt.show()
  
  
  # In[11]:
  
  
  do_cross_validation(X,y)
  
  
  # In[12]:
  
  
  from keras.models import load_model
  X,y = prepare_inputs_len_x_alphabet(pc_test,nc_test,ALPHABET)
  best_model=load_model(MODELPATH)
  scores = best_model.evaluate(X, y, verbose=0)
  print("The best model parameters were saved during cross-validation.")
  print("Best was defined as maximum validation accuracy at end of any epoch.")
  print("Now re-load the best model and test it on previously unseen data.")
  print("Test on",len(pc_test),"PC seqs")
  print("Test on",len(nc_test),"NC seqs")
  print("%s: %.2f%%" % (best_model.metrics_names[1], scores[1]*100))
  
  
  # In[13]:
  
  
  from sklearn.metrics import roc_curve
  from sklearn.metrics import roc_auc_score
  ns_probs = [0 for _ in range(len(y))]
  bm_probs = best_model.predict(X)
  ns_auc = roc_auc_score(y, ns_probs)
  bm_auc = roc_auc_score(y, bm_probs)
  ns_fpr, ns_tpr, _ = roc_curve(y, ns_probs)
  bm_fpr, bm_tpr, _ = roc_curve(y, bm_probs)
  plt.plot(ns_fpr, ns_tpr, linestyle='--', label='Guess, auc=%.4f'%ns_auc)
  plt.plot(bm_fpr, bm_tpr, marker='.', label='Model, auc=%.4f'%bm_auc)
  plt.title('ROC')
  plt.xlabel('False Positive Rate')
  plt.ylabel('True Positive Rate')
  plt.legend()
  plt.show()
  print("%s: %.2f%%" %('AUC',bm_auc*100.0))
  
  
  # In[14]:
  
  
  t = time.time()
  time.strftime('%Y-%m-%d %H:%M:%S %Z', time.localtime(t))
  
  
  # In[14]:
  
  
  
  
  