set half precision for numpy arrays

scripts/genome_wide/add_win_channels.py

@@ -17,7 +17,7 @@ def init_log(logfile):
 
 
 def parse_args():
-    default_win = 25
+    default_win = 200
     parser = argparse.ArgumentParser(
         description='Add window specific channels')
 
@@ -187,7 +187,7 @@ def get_reads(chrom, pos):
     # get starting time
     last_t = time()
     # Initialize numpy array
-    X_enh = np.zeros(shape=(X.shape[:2] + (len(ch),)), dtype=np.int8)
+    X_enh = np.zeros(shape=(X.shape[:2] + (len(ch),)), dtype=np.half)
 
     for i, p in enumerate(y.keys(), start=0):
         # Every n_r alignments, write log informations
@@ -245,6 +245,7 @@ def main():
 
     with pysam.AlignmentFile(args.bam, "rb") as bam:
         X, y = add_channels(args, bam)
+        X = X.astype(np.half)
         save_windows(X, y, args.output)
     logging.info('Finished in %f seconds' % (time() - t0))
 

scripts/genome_wide/chr_array.py

@@ -93,7 +93,7 @@ def create_carray(ibam, chrom, twobit, bigwig, outDir, cmd_name):
     chrlen = get_chr_len(ibam, chrom)
     channel_index = 0
     n_channels = 54
-    chr_array = np.zeros(shape=(chrlen, n_channels), dtype=np.float64)
+    chr_array = np.zeros(shape=(chrlen, n_channels), dtype=np.half)
     bw_map = pyBigWig.open(bigwig)
     # dictionary of key choices
     direction_list = {
@@ -147,7 +147,7 @@ def create_carray(ibam, chrom, twobit, bigwig, outDir, cmd_name):
                                       dtype=int)
                     vals = np.fromiter(channel_data[chrom][current_channel]
                                        [split_direction].values(),
-                                       dtype=np.float32)
+                                       dtype=np.half)
                     if len(idx) > 0:
                         chr_array[idx, channel_index] = vals
                         assert chr_array[idx, channel_index].any(), \
@@ -205,7 +205,7 @@ def create_carray(ibam, chrom, twobit, bigwig, outDir, cmd_name):
                  (current_channel, channel_index))
     chr_array[:, channel_index] = np.array(bw_map.values(
         chrom, 0, chrlen),
-        dtype=np.float32)
+        dtype=np.half)
     channel_index += 1
     current_channel = 'one_hot_encoding'
     logging.info("Adding channel %s at index %d" %
@@ -222,7 +222,7 @@ def create_carray(ibam, chrom, twobit, bigwig, outDir, cmd_name):
 
 
 def main():
-    default_chr = '22'
+    default_chr = '12'
     parser = argparse.ArgumentParser(
         description='Create channels from saved data')
     parser.add_argument('-b',

scripts/genome_wide/clipped_read_distance.py

@@ -102,7 +102,7 @@ def main():
     parser.add_argument('-c',
                         '--chr',
                         type=str,
-                        default='12',
+                        default='22',
                         help="Specify chromosome")
     parser.add_argument('-o',
                         '--out',

scripts/genome_wide/coverage.py

@@ -46,7 +46,7 @@ def get_coverage(ibam, chrName, minMAPQ, outFile):
     chrLen = [i['LN'] for i in header_dict['SQ'] if i['SN'] == chrName][0]
     start_pos = 0
     stop_pos = chrLen
-    cov = np.zeros((chrLen, 5))
+    cov = np.zeros((chrLen, 5), dtype=np.int)
     n_r = 10 ** 6
     last_t = time()
 

scripts/genome_wide/create_window_pairs.py

@@ -23,6 +23,7 @@ def load_chr_array(channel_data_dir, chrlist):
         logging.info("Loading file %s" % carray_file)
         chr_array[c] = bcolz.open(rootdir=carray_file)
         logging.info("Array shape: %s" % str(chr_array[c].shape))
+        logging.info("Array data type: %s" % str(chr_array[c].dtype))
     return chr_array
 
 
@@ -74,7 +75,7 @@ def get_windows(carrays_dir, outDir, chrom_list, win, label_file_path, mode, npz
         n_r = 10 ** 5
         last_t = time()
         i = 1
-        padding = np.zeros(shape=(padding_len, n_channels), dtype=np.float32)
+        padding = np.zeros(shape=(padding_len, n_channels), dtype=np.half)
         if npz_mode:
             numpy_array = []
         logging.info('Creating np.arrays win1 and win2...')
@@ -105,6 +106,7 @@ def get_windows(carrays_dir, outDir, chrom_list, win, label_file_path, mode, npz
             for i in np.arange(numpy_array.shape[2]):
                 logging.info("windows array: non-zero elements at index %d:%d" %
                              (i, np.argwhere(numpy_array[i, :] != 0).shape[0]))
+            numpy_array= numpy_array.astype(np.half)
             np.savez(file=os.path.join(outDir, 'windows'),
                      data=numpy_array,
                      labels=labs)
@@ -131,7 +133,7 @@ def main():
     parser.add_argument('-p',
                         '--outputpath',
                         type=str,
-                        default='./cnn/win25/split_reads/windows/DEL',
+                        default='./cnn/win200/split_reads/windows/DEL',
                         help="Specify output path")
     parser.add_argument('-l',
                         '--logfile',
@@ -140,12 +142,12 @@ def main():
     parser.add_argument('-w',
                         '--window',
                         type=int,
-                        default=25,
+                        default=200,
                         help="Specify window size")
     parser.add_argument('-lb',
                         '--labels',
                         type=str,
-                        default='./cnn/win25/split_reads/windows/DEL/labels.json.gz',
+                        default='./cnn/win200/split_reads/windows/DEL/labels.json.gz',
                         help="Specify label file")
     parser.add_argument('-m',
                         '--mode',

scripts/genome_wide/model_functions.py

@@ -17,63 +17,7 @@ def unfold_win_id(win_id):
     return chr1, pos1, chr2, pos2, strand_info
 
 
-# def create_model_with_mcfly(X, y_binary):
-#
-#     models = modelgen.generate_models(X.shape,
-#                                       y_binary.shape[1],
-#                                       number_of_models=1,
-#                                       model_type='CNN',
-#                                       cnn_min_layers=2,
-#                                       cnn_max_layers=2,
-#                                       cnn_min_filters=4,
-#                                       cnn_max_filters=4,
-#                                       cnn_min_fc_nodes=6,
-#                                       cnn_max_fc_nodes=6,
-#                                       low_lr=4, high_lr=4,
-#                                       low_reg=1, high_reg=1,
-#                                       kernel_size=7)
-#
-#     i = 0
-#     for model, params, model_types in models:
-#         logging.info('model ' + str(i))
-#         i = i + 1
-#         logging.info(params)
-#         logging.info(model.summary())
-#
-#     return models
-#
-#
-# def train_model_with_mcfly(model, xtrain, ytrain, xval, yval):
-#
-#     train_set_size = xtrain.shape[0]
-#     nr_epochs = 5
-#
-#     histories, val_accuracies, val_losses = find_architecture.train_models_on_samples(xtrain, ytrain,
-#                                                                                       xval, yval,
-#                                                                                       model, nr_epochs=nr_epochs,
-#                                                                                       subset_size=train_set_size,
-#                                                                                       verbose=True)
-#
-#     best_model_index = np.argmax(val_accuracies)
-#     best_model, best_params, best_model_types = model[best_model_index]
-#     # print(best_model_index, best_model_types, best_params)
-#
-#     history = best_model.fit(xtrain, ytrain,
-#                              epochs=nr_epochs,
-#                              validation_data=(xval, yval),
-#                              verbose=True,
-#                              shuffle=True)
-#
-#     return history, best_model
-
-def get_data(windows_list, npz_mode, svtype):
-
-    def filter_labels(X, y, win_ids):
-        keep = [i for i, v in enumerate(y) if v in [svtype, 'no' + svtype]]
-        X = X[np.array(keep)]
-        y = [y[i] for i in keep]
-        win_ids = [win_ids[i] for i in keep]
-        return X, y, win_ids
+def get_data(windows_list, svtype):
 
     X = []
     y = []

scripts/genome_wide/snv.py

@@ -32,7 +32,7 @@ def get_snv_number(query_seq_list, reference_base):
     stop_pos = chrLen
     reference_sequence = twobit.TwoBitFile(itwobit)
     snv_list = ['BQ', 'SNV', 'MAPQ']
-    snv_array = np.zeros(shape=(stop_pos, len(snv_list)), dtype=np.float32)
+    snv_array = np.zeros(shape=(stop_pos, len(snv_list)), dtype=np.half)
     snv_dict = {v: n for n, v in enumerate(snv_list)}
 
     for pileupcolumn in bamfile.pileup(chrName,

scripts/genome_wide/train.py

@@ -83,11 +83,11 @@ def create_model(dim_length, dim_channels, outputdim):
     model.add(Activation('relu'))  # Relu activation
 
     # Adding one more FC layer
-    model.add(
-        Dense(units=model_params['fc_nodes'],
-              kernel_regularizer=l2(regularization_rate),
-              kernel_initializer=weightinit))  # Fully connected layer
-    model.add(Activation('relu'))  # Relu activation
+    # model.add(
+    #     Dense(units=model_params['fc_nodes'],
+    #           kernel_regularizer=l2(regularization_rate),
+    #           kernel_initializer=weightinit))  # Fully connected layer
+    # model.add(Activation('relu'))  # Relu activation
 
     model.add(Dense(units=outputdim, kernel_initializer=weightinit))
     model.add(BatchNormalization())
@@ -156,10 +156,10 @@ def train(model_fn, params, X_train, y_train, y_train_binary):
         epochs=model_params['epochs'],
         shuffle=True,
         class_weight=class_weights,
-        verbose=1,
+        verbose=0,
         callbacks=callbacks)
 
-    model = load_model(model_fn)
+    # model = load_model(model_fn)
 
     return model, history, X_train.shape[0], int(X_train.shape[0] *
                                                  model_params['validation_split'])
@@ -199,8 +199,9 @@ def cv_train_and_evaluate(X, y, y_binary, win_ids, train_indices, test_indices,
     results.to_csv(os.path.join(model_dir, 'metrics.csv'), sep='\t')
 
 
-def cross_validation(training_windows, outDir, npz_mode, svtype, kfold):
-    X, y, win_ids = get_data(training_windows, npz_mode, svtype)
+def cross_validation(training_windows, outDir, svtype, kfold):
+
+    X, y, win_ids = get_data(training_windows, svtype)
     y_binary = to_categorical(y, num_classes=len(mapclasses.keys()))
 
     # Instantiate the cross validator
@@ -217,8 +218,8 @@ def cross_validation(training_windows, outDir, npz_mode, svtype, kfold):
                               train_indices, test_indices, model_dir, svtype)
 
 
-def cross_validation_by_chrom(training_windows, outDir, npz_mode, svtype, chrlist):
-    X, y, win_ids = get_data(training_windows, npz_mode, svtype)
+def cross_validation_by_chrom(training_windows, outDir, svtype, chrlist):
+    X, y, win_ids = get_data(training_windows, svtype)
     y_binary = to_categorical(y, num_classes=len(mapclasses.keys()))
 
     # print(win_ids)
@@ -250,10 +251,10 @@ def cross_validation_by_chrom(training_windows, outDir, npz_mode, svtype, chrlis
 
 
 def train_and_test_model(training_name, test_name, training_windows, test_windows,
-                         outDir,
-                         npz_mode, svtype):
-    X_test, y_test, win_ids_test = get_data(training_windows, npz_mode, svtype)
-    X_test, y_test, win_ids_test = get_data(test_windows, npz_mode, svtype)
+                         outDir, svtype):
+
+    X_train, y_train, win_ids_train = get_data(training_windows, svtype)
+    X_test, y_test, win_ids_test = get_data(test_windows, svtype)
 
     # Parameters
     params = {
@@ -346,12 +347,12 @@ def main():
     parser.add_argument('-e',
                         '--epochs',
                         type=int,
-                        default=1,
+                        default=10,
                         help="Number of epochs")
     parser.add_argument('-b',
                         '--batch_size',
                         type=int,
-                        default=32,
+                        default=4,
                         help="Batch size")
     parser.add_argument('-val',
                         '--validation_split',
@@ -366,22 +367,22 @@ def main():
     parser.add_argument('-cnn_layers',
                         '--cnn_layers',
                         type=int,
-                        default=4,
+                        default=2,
                         help="Number of convolutional layers")
     parser.add_argument('-cnn_filters',
                         '--cnn_filters',
                         type=int,
-                        default=8,
+                        default=32,
                         help="Number of convolutional filters")
     parser.add_argument('-kernel_size',
                         '--kernel_size',
                         type=int,
-                        default=7,
+                        default=5,
                         help="Number of convolutional filters")
     parser.add_argument('-fc_nodes',
                         '--fc_nodes',
                         type=int,
-                        default=16,
+                        default=6,
                         help="Number of neurons in the dense layer")
     parser.add_argument('-learning_rate_exp',
                         '--learning_rate_exp',
@@ -433,22 +434,19 @@ def main():
             training_name=args.training_sample_name,
             test_name=args.test_sample_name,
             training_windows=training_windows_list,
-            test_windows=args.test_windowsr,
+            test_windows=args.test_windows,
             outDir=output_dir,
-            npz_mode=args.load_npz,
             svtype=args.svtype
         )
     else:
         if args.cv == 'kfold':
             cross_validation(training_windows=training_windows_list,
                              outDir=output_dir,
-                             npz_mode=args.load_npz,
                              svtype=args.svtype,
                              kfold=args.kfold)
         elif args.cv == 'chrom':
             cross_validation_by_chrom(training_windows=training_windows_list,
                                       outDir=output_dir,
-                                      npz_mode=args.load_npz,
                                       svtype=args.svtype,
                                       chrlist=args.chrlist)
     logging.info('Elapsed time training and testing = %f seconds' %