<reponame>genepii/seqmet
from __future__ import print_function
import os
import sys
import getopt
import math
import random
import numpy as np

##Generate oligonucleotides with a maximal hamming distance
#v0.0.1

def main(argv):
    global oup
    global gc
    global length
    global size
    oup = 'genoligos.txt'
    gc = [35,65]
    length = 150
    size = 800
    try:
        opts, args = getopt.getopt(argv, 'hi:o:g:l:s:', ['help', 'output', 'gc', 'length', 'size'])
        for opt, arg in opts:
            if opt in ('-h', '--help'):
                usage()
                sys.exit()
            elif opt in ('-o', '--output'):
                oup = arg
            elif opt in ('-g', '--gc'):
                gc = [int(arg.split(',')[0]), int(arg.split(',')[1])]
            elif opt in ('-l', '--length'):
                length = int(arg)
            elif opt in ('-s', '--size'):
                size = int(arg)
    except getopt.GetoptError:
        usage()
        sys.exit(2)

def usage():
    print('usage: ' + sys.argv[0] + ' -h --help -o --output [txt] -g --gc [integer] -l --length [integer] -s --size [integer]')

if __name__ == '__main__':
    main(sys.argv[1:])

max_gc = int(math.floor((float(gc[1])/100)*length))
min_gc = int(math.ceil((float(gc[0])/100)*length))

base = ['A', 'T', 'G', 'C']
oligos = []
temp_oligos = ['A', 'T', 'G', 'C']

def averageOfList(numOfList):
    avg = sum(numOfList) / len(numOfList)
    return avg

def add_rd_nt(ol, nt):
    if len(ol) < 2:
        return ol + nt
    if ((ol[len(ol)-1] != nt or ol[len(ol)-2] != nt)) and (((ol + nt).count('G') + (ol + nt).count('C'))) <= max_gc and (ol + nt).count('ATATAT') == 0 and (ol + nt).count('TATATA') == 0 and (ol + nt).count('ACACAC') == 0 and (ol + nt).count('CACACA') == 0 and (ol + nt).count('AGAGAG') == 0 and (ol + nt).count('GAGAGA') == 0 and (ol + nt).count('TGTGTG') == 0 and (ol + nt).count('GTGTGT') == 0 and (ol + nt).count('TCTCTC') == 0 and (ol + nt).count('CTCTCT') == 0 and (ol + nt).count('GCGCGC') == 0 and (ol + nt).count('CGCGCG') == 0:
        return ol + nt
    else:
        return False

def hamming_distance(s1, s2):
    #Return the Hamming distance between equal-length sequences
    if len(s1) != len(s2):
        raise ValueError("Undefined for sequences of unequal length")
    return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))

for i in range(length):

    sys.stdout.flush()
    print("Selecting oligos of " + str(i+1) + "nt length among " + str(len(temp_oligos)) + " sequences", end='\r')

    oligos = []

    if len(temp_oligos) < size:
        current_size = len(temp_oligos)
    else:
        current_size = size

    while len(oligos) < current_size:
        hamming_oligos = []
        for j in range(len(temp_oligos)):
            temp_hamming_oligos = []
            if len(oligos) >= 1:
                for h in range(len(oligos)):
                    temp_hamming_oligos.append(hamming_distance(temp_oligos[j], oligos[h]))
                hamming_oligos.append(min(temp_hamming_oligos))
            else:
                hamming_oligos.append(100)

        best_idx = np.argpartition(np.array(hamming_oligos), -(1))[-(1):]
        oligos.append(temp_oligos.pop(best_idx))

    temp_oligos = []
    for j in range(len(oligos)):
        random.shuffle(base)
        if current_size < size:
            for k in range(len(base)):
                new_oligos = add_rd_nt(oligos[j], base[k])
                if new_oligos is not False:
                    temp_oligos.append(new_oligos)
        else:
            for k in range(2):
                new_oligos = add_rd_nt(oligos[j], base[k])
                if new_oligos is not False:
                    temp_oligos.append(new_oligos)

    #if len(temp_oligos) > size*2:
    #    random.shuffle(temp_oligos)
    #    temp_oligos = temp_oligos[:size*2]


oligos_final = []

for i in range(len(oligos)):
    if (oligos[i].count('G') + oligos[i].count('C')) >= 0:
        oligos_final.append(oligos[i])

oligos_minhd = []
oligos_maxhd = []
oligos_meanhd = []
for i in range(len(oligos_final)):
    temp_hamming_oligos = []
    for j in range(len(oligos_final)):
        if i != j:
            temp_hamming_oligos.append(hamming_distance(oligos_final[i], oligos_final[j]))
    oligos_minhd.append(min(temp_hamming_oligos))
    oligos_maxhd.append(max(temp_hamming_oligos))
    oligos_meanhd.append(averageOfList(temp_hamming_oligos))

print(oligos_final)

w = open(oup, 'w')
for i in range(len(oligos_final)):
    w.write('>oligo_' + str(i) + '|min_hd=' + str(oligos_minhd[i]) + '\n' + oligos_final[i] + '\n')
w.close()

print('total_mean_hd=' + str(averageOfList(oligos_meanhd)) + ' total_min_hd=' + str(min(oligos_minhd)) + ' total_max_hd=' + str(max(oligos_maxhd)) + '                   \n')

from __future__ import print_function
import os
import sys
import getopt
import math
import random
import numpy as np

##Generate oligonucleotides with a maximal hamming distance
#v0.0.1

def main(argv):
    global oup
    global gc
    global length
    global size
    oup = 'genoligos.txt'
    gc = [35,65]
    length = 150
    size = 800
    try:
        opts, args = getopt.getopt(argv, 'hi:o:g:l:s:', ['help', 'output', 'gc', 'length', 'size'])
        for opt, arg in opts:
            if opt in ('-h', '--help'):
                usage()
                sys.exit()
            elif opt in ('-o', '--output'):
                oup = arg
            elif opt in ('-g', '--gc'):
                gc = [int(arg.split(',')[0]), int(arg.split(',')[1])]
            elif opt in ('-l', '--length'):
                length = int(arg)
            elif opt in ('-s', '--size'):
                size = int(arg)
    except getopt.GetoptError:
        usage()
        sys.exit(2)

def usage():
    print('usage: ' + sys.argv[0] + ' -h --help -o --output [txt] -g --gc [integer] -l --length [integer] -s --size [integer]')

if __name__ == '__main__':
    main(sys.argv[1:])

max_gc = int(math.floor((float(gc[1])/100)*length))
min_gc = int(math.ceil((float(gc[0])/100)*length))

base = ['A', 'T', 'G', 'C']
oligos = []
temp_oligos = ['A', 'T', 'G', 'C']

def averageOfList(numOfList):
    avg = sum(numOfList) / len(numOfList)
    return avg

def add_rd_nt(ol, nt):
    if len(ol) < 2:
        return ol + nt
    if ((ol[len(ol)-1] != nt or ol[len(ol)-2] != nt)) and (((ol + nt).count('G') + (ol + nt).count('C'))) <= max_gc and (ol + nt).count('ATATAT') == 0 and (ol + nt).count('TATATA') == 0 and (ol + nt).count('ACACAC') == 0 and (ol + nt).count('CACACA') == 0 and (ol + nt).count('AGAGAG') == 0 and (ol + nt).count('GAGAGA') == 0 and (ol + nt).count('TGTGTG') == 0 and (ol + nt).count('GTGTGT') == 0 and (ol + nt).count('TCTCTC') == 0 and (ol + nt).count('CTCTCT') == 0 and (ol + nt).count('GCGCGC') == 0 and (ol + nt).count('CGCGCG') == 0:
        return ol + nt
    else:
        return False

def hamming_distance(s1, s2):
    #Return the Hamming distance between equal-length sequences
    if len(s1) != len(s2):
        raise ValueError("Undefined for sequences of unequal length")
    return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))

for i in range(length):

    sys.stdout.flush()
    print("Selecting oligos of " + str(i+1) + "nt length among " + str(len(temp_oligos)) + " sequences", end='\r')

    oligos = []

    if len(temp_oligos) < size:
        current_size = len(temp_oligos)
    else:
        current_size = size

    while len(oligos) < current_size:
        hamming_oligos = []
        for j in range(len(temp_oligos)):
            temp_hamming_oligos = []
            if len(oligos) >= 1:
                for h in range(len(oligos)):
                    temp_hamming_oligos.append(hamming_distance(temp_oligos[j], oligos[h]))
                hamming_oligos.append(min(temp_hamming_oligos))
            else:
                hamming_oligos.append(100)

        best_idx = np.argpartition(np.array(hamming_oligos), -(1))[-(1):]
        oligos.append(temp_oligos.pop(best_idx))

    temp_oligos = []
    for j in range(len(oligos)):
        random.shuffle(base)
        if current_size < size:
            for k in range(len(base)):
                new_oligos = add_rd_nt(oligos[j], base[k])
                if new_oligos is not False:
                    temp_oligos.append(new_oligos)
        else:
            for k in range(2):
                new_oligos = add_rd_nt(oligos[j], base[k])
                if new_oligos is not False:
                    temp_oligos.append(new_oligos)

    #if len(temp_oligos) > size*2:
    #    random.shuffle(temp_oligos)
    #    temp_oligos = temp_oligos[:size*2]


oligos_final = []

for i in range(len(oligos)):
    if (oligos[i].count('G') + oligos[i].count('C')) >= 0:
        oligos_final.append(oligos[i])

oligos_minhd = []
oligos_maxhd = []
oligos_meanhd = []
for i in range(len(oligos_final)):
    temp_hamming_oligos = []
    for j in range(len(oligos_final)):
        if i != j:
            temp_hamming_oligos.append(hamming_distance(oligos_final[i], oligos_final[j]))
    oligos_minhd.append(min(temp_hamming_oligos))
    oligos_maxhd.append(max(temp_hamming_oligos))
    oligos_meanhd.append(averageOfList(temp_hamming_oligos))

print(oligos_final)

w = open(oup, 'w')
for i in range(len(oligos_final)):
    w.write('>oligo_' + str(i) + '|min_hd=' + str(oligos_minhd[i]) + '\n' + oligos_final[i] + '\n')
w.close()

print('total_mean_hd=' + str(averageOfList(oligos_meanhd)) + ' total_min_hd=' + str(min(oligos_minhd)) + ' total_max_hd=' + str(max(oligos_maxhd)) + '                   \n')
##Generate oligonucleotides with a maximal hamming distance #v0.0.1 #Return the Hamming distance between equal-length sequences #if len(temp_oligos) > size*2: #    random.shuffle(temp_oligos) #    temp_oligos = temp_oligos[:size*2]