IGEM:Harvard/2006/Container Design 4/Python Code/Full Code Final

From OpenWetWare

Jump to: navigation, search

Modified Versions of William's Code - Finalized 7/11/06

Main

#!/usr/bin/python

import sys
# for regular expression matching
import re
# for reading structures in from files
import pickle
from honeycomb_pointers_v1 import *
from pointers_to_sequences_v1 import *

filename = raw_input('Enter the filename for the text-based node lattice array: ')
num_zones = int(raw_input('Enter the number of 42bp zones per double helix: '))
#filename = 'lid_1_ascii.txt'
#num_zones = 3
periodic_structure_flag = False

node_ra = read_text_format_node_array(filename)
print_node_lattice_array(node_ra)
TPP_ra = token_pointer_pair_array(node_ra, num_zones, periodic_structure_flag)
path_ra = token_pointer_path_array(TPP_ra)
OTP_ra = oligo_token_pointer_array(path_ra)
check_token_representation(OTP_ra, TPP_ra)


# Print out the longest path, just for fun
longest_path = []
for path in path_ra:
	if len(path) > len(longest_path):
		longest_path = path
print_path(longest_path, TPP_ra)

#####
# Oligo splitting - this time reading from a file and not asking for user
# input
#####
fin_barrel = None
fin_lid = None

try:
	fin_barrel = open("barrel_oligos_to_split.txt", "r")
	fin_lid = open("lid_oligos_to_split.txt", "r")
except IOError, e:
	print "Error in file IO: ", e

# Ask the user if they are running a lid or a barrel

shape = int(raw_input("Enter 1 if you are running a barrel, 2 if lid: "))
if (shape == 1):
	oligos_to_split = pickle.load(fin_barrel)
elif (shape == 2):
	oligos_to_split = pickle.load(fin_lid)
else:
	print 'Please modify code or run with lid or 30hb barrel'
	

new_OTP_ra = OTP_ra[:]
for pair in oligos_to_split:
	oligo_num = pair[0]
	print oligo_num
	num_toks = pair[1]
	print num_toks
	
       	new_OTP_ra = split_oligo(new_OTP_ra, oligo_num, num_toks)

       	print new_OTP_ra
     	print len(OTP_ra)
       	print len(new_OTP_ra)

# if it's the barrel design and so that all the numbers aren't messed
# up, the 7bp token on the start of strand 10 needs to be removed
# because it's going to be left unpaired

if (shape == 1):
	new_OTP_ra = new_OTP_ra[:61] + new_OTP_ra[62:]

if fin_barrel: fin_barrel.close()
if fin_lid: fin_lid.close()




####
# generate and print the oligo token grid
####

# Initialize the grid with all periods
num_strands = len(TPP_ra)
num_subzones = len(TPP_ra[0])

sub_token_visit_ra = ['.' for subzone_num in range(num_subzones)]
grid_ra = [sub_token_visit_ra[:] for strand_num in range(num_strands)]

oligo_num = 0
for oligo in new_OTP_ra:
	grid_ra = generate_oligo_path(oligo, oligo_num, grid_ra) 
	oligo_num = oligo_num + 1
print grid_ra

print_all_oligos(grid_ra, num_strands, num_subzones)

####
# Generate oligo sequences
####
strand_ra = strand_array('lid_1_scaffold.txt', TPP_ra)
token_ra = token_array(strand_ra)
oligo_ra = oligo_array(new_OTP_ra, token_ra)

# The array of oligos that does have the oligos split but does not have the
# aptamers or latches added.
original_oligo_ra = oligo_ra[:]
latch_2_oligo_ra = oligo_ra[:]

#####
# Add apts this time using file input instead of user input
#####

# Constants
apt_seq = 'GGTTGGTGTGGTTGG'
T_linker = 'TTT'

fin_barrel = None
        
try:
        fin_barrel = open("barrel_apts_to_add.txt", "r")
except IOError, e:
        print "Error in file IO: ", e
        
if (shape == 1):
    	apts_to_add = pickle.load(fin_barrel)

	for apt_specs in apts_to_add:
		oligo_num = apt_specs[0]
		type = apt_specs[1]
		if (type == 1):
			# apt is pointing in so add 'I' as a flag at the end
			oligo_ra[oligo_num] = oligo_ra[oligo_num] + T_linker + apt_seq + 'I'
		elif (type == 2):
			# apt is pointing out so add 'O' as a flag
			oligo_ra[oligo_num] = oligo_ra[oligo_num] + T_linker + apt_seq + 'O'
		else:
			# incorrect type
			print 'Bad input - aptamer needs to be pointing in or out'

#####
# Add latches
#####

fin_barrel = None
fin_lid1 = None
fin_lid2 = None
fin_barrel_design2 = None
fin_lid1_design2 = None
fin_lid2_design2 = None
                
try:
        fin_barrel = open("barrel_latches_to_add.txt", "r")
	fin_lid1 = open("lid1_latches_to_add.txt", "r")
	fin_lid2 = open("lid2_latches_to_add.txt", "r")
	fin_barrel_design2 = open("barrel_latches_to_add_design2.txt", "r")
	fin_lid1_design2 = open("lid1_latches_to_add_design2.txt", "r")
	fin_lid2_design2 = open("lid2_latches_to_add_design2.txt", "r")
except IOError, e:
        print "Error in file IO: ", e


if (shape == 1):
        latches_to_add = pickle.load(fin_barrel)
	latches_to_add_d2 = pickle.load(fin_barrel_design2)
elif (shape == 2):
	# Ask which lid is being run
	print 'What lid are you running?'
	type = int(raw_input('Enter 1 if lid1, 2 if lid2: '))
	if (type == 1):
		latches_to_add = pickle.load(fin_lid1)
		latches_to_add_d2 = pickle.load(fin_lid1_design2)
	elif (type == 2):
		latches_to_add = pickle.load(fin_lid2)
		latches_to_add_d2 = pickle.load(fin_lid2_design2)

        
# note: latch sequence includes any linker Ts already
for latch_specs in latches_to_add:
	oligo_num = latch_specs[0]
        latch = latch_specs[1]
	oligo_ra[oligo_num] = oligo_ra[oligo_num] + latch

for latch_specs in latches_to_add_d2:
	oligo_num = latch_specs[0]
	latch = latch_specs[1]
	latch_2_oligo_ra[oligo_num] = latch_2_oligo_ra[oligo_num] + latch


#####
# oligo sorting
#####
#
barrel_core = [0, 3, 4, 5, 6, 8, 9, 10, 12, 13, 14, 15, 16, 18, 19, 20, 22, 24, 25, 26, 28, 29, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77]
barrel_latch = [21, 27, 43, 60] 
barrel_aptamer_out = [2, 7, 17, 23] 
barrel_aptamer_in = [1, 11, 36, 47]
lid_core = [0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 30, 31, 32, 33, 34]
lid_latch = [2, 29]

# output all oligos
output_file = file('lid_1_all_oligos.txt', 'w')
for oligo in oligo_ra:
	output_file.write(oligo + '\n')
output_file.close()

if (shape == 1):
	# barrel output
	
	barrel_oligos_file = file('barrel_oligos_sorted.txt', 'w')
	# output barrel core oligos
	barrel_oligos_file.write('Core Oligos' + '\n')
	barrel_oligos_file.close()
	# to append the file needs to be reopened in append mode.
	barrel_oligos_file = file('barrel_oligos_sorted.txt', 'a')

	for oligo_num in barrel_core:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')

	# out + apts
	barrel_oligos_file.write('\n' + 'Out + Apts' + '\n')
	for oligo_num in barrel_aptamer_out:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')
	# out - apts
	barrel_oligos_file.write('\n' + 'Out - Apts' + '\n')
	for oligo_num in barrel_aptamer_out:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + original_oligo_ra[oligo_num] + '\n')
	
	# in - apts
	barrel_oligos_file.write('\n' + 'In - Apts' + '\n')
	for oligo_num in barrel_aptamer_in:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + original_oligo_ra[oligo_num] + '\n')

	# in + apts
	barrel_oligos_file.write('\n' + 'In + Apts' +'\n')
	for oligo_num in barrel_aptamer_in:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')
	
	# barrel oligos with latches - two oligo design
	barrel_oligos_file.write('\n' + 'Oligos with Latches - Design 1' + '\n')
	for oligo_num in barrel_latch:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')
	
	# barrel oligos with latches - four oligo design
	barrel_oligos_file.write('\n' + 'Oligos with Latches - Design 2\n')
	for oligo_num in barrel_latch:
		barrel_oligos_file.write(str(oligo_num) + ' Design 2: ' + latch_2_oligo_ra[oligo_num] + '\n')

	# barrel oligos replacing all latches
	barrel_oligos_file.write('\n' + 'Oligos replacing all latches\n')
	for oligo_num in barrel_latch:
		barrel_oligos_file.write(str(oligo_num) + ' : ' + original_oligo_ra[oligo_num] + '\n')

	
	barrel_oligos_file.close()


if (shape == 2):
	# lid output

	# output lid core oligos
	# we first open the file for writing because we want to clear 
	# anything that was there before. Then we close it and reopen it
	# in append mode.
	lid_oligos_file = file('lid1_oligos_sorted.txt', 'w')
	lid_oligos_file.write('Lid Core' + '\n')
	lid_oligos_file.close()
	lid_oligos_file = file('lid1_oligos_sorted.txt', 'a')
	for oligo_num in lid_core:
		lid_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')

	# output lid oligos with latches -- 2 oligo design
	lid_oligos_file.write('\nLid Latch Oligos - 2 Oligo Design\n')
	for oligo_num in lid_latch:
		lid_oligos_file.write(str(oligo_num) + ' : ' + oligo_ra[oligo_num] + '\n')

	# output lid oligos with latches -- 4 oligo design
	lid_oligos_file.write('\nLid Latch Oligos - 4 Oligo Design\n')
	for oligo_num in lid_latch:
		lid_oligos_file.write(str(oligo_num) + ' Design 2 : ' + latch_2_oligo_ra[oligo_num] + '\n')

	# output lid oligos that replace those with latches
	lid_oligos_file.write('\nLid Replacement Oligos - No Latches\n')
	for oligo_num in lid_latch:
		lid_oligos_file.write(str(oligo_num) + ' : ' + original_oligo_ra[oligo_num] + '\n')
	
	lid_oligos_file.close()	

sys.stdout.write('\n\n\n\n\n\n')

Honeycomb_pointers_v1.py

#!/usr/bin/python

import sys






######
# This function reads in the node array from a text file
# It adds a border of nodes automatically
######
def read_text_format_node_array(filename):
	# Read in file
	input_file = file(filename, 'r')
	lines = input_file.readlines()
	input_file.close()
	row_string_ra = [line[:-1] for line in lines]

	
	# Check to make sure each line is the same length
	no_length_violation = True
	length = len(row_string_ra[0])
	for row_string in row_string_ra:
		if len(row_string) != length and len(row_string) != 0:
			sys.stdout.write('ERROR: Not all lines of inputted node lattice array are the same length.\n')
			sys.stdout.write('Length is ' + str(len(row_string)) + '\n')
			no_length_violation = False

		
	# Parse into pre node array
	pre_node_ra = []
	for row_string in row_string_ra:
		num_row_nodes = len(row_string)/3
		sub_pre_node_ra = []
		for row_node_num in range(num_row_nodes):
			sub_pre_node_ra.append(row_string[row_node_num*3:row_node_num*3 + 3])
		if sub_pre_node_ra != []:
			pre_node_ra.append(sub_pre_node_ra)	

	
	# Parse pre node array into node array
	num_rows = len(pre_node_ra)/2
	num_row_nodes = len(pre_node_ra[0])
	node_ra = [['.' for row_node_num in range(num_row_nodes + num_row_nodes%2 + 2)]]
	for row_num in range(num_rows):
		sub_node_ra = ['.']
		for row_node_num in range(num_row_nodes):
			pre_node_string = pre_node_ra[row_num*2 + (row_node_num + row_num)%2][row_node_num]
			if pre_node_string == '...':
				node = '.'
			else:
				node = int(pre_node_string)
			sub_node_ra.append(node)
		sub_node_ra.append('.')
		if num_row_nodes%2 == 1:
			sub_node_ra.append('.')			
		node_ra.append(sub_node_ra)
	node_ra.append(['.' for row_node_num in range(num_row_nodes + num_row_nodes%2 + 2)])

	
	# Check for parity violations
	num_parity_violations = 0
	for row_num in range(num_rows):
		for row_node_num in range(num_row_nodes):
			node = node_ra[row_num][row_node_num]
			if node != '.':
				if (node + row_num + row_node_num)%2 == 1:
					sys.stdout.write('ERROR: Parity violation for strand ' + str(node) + '.\n')
					sys.stdout.write('Parity is the row number plus the row-node number.\n')
					sys.stdout.write('Make sure even-numbered strands are on even-parity nodes.\n')
					sys.stdout.write('Make sure odd-numbered strands are on odd-parity nodes.\n')
					num_parity_violations += 1

	
	if num_parity_violations > 0 and no_length_violations == True:
		return
	else:
		return node_ra






######
# This function prints the node lattice array in a honeycomb format
######
def print_node_lattice_array(node_ra):
	sys.stdout.write('\nNODE LATTICE ARRAY\n')
	zeroes = '000'
	for y in range(len(node_ra)):
		even_row_string = ''
		odd_row_string = '   '
		for x in range(len(node_ra[0])):
			string_element = str(node_ra[y][x])
			if string_element == '.':
				string_element += '..'
			else:
				string_element = zeroes[:3 - len(string_element)] + string_element
			if x%2 == 0:
				even_row_string += string_element + '   '
			else:
				odd_row_string += string_element + '   '
		if y%2 == 0:
			sys.stdout.write(even_row_string + '\n')
			sys.stdout.write(odd_row_string + '\n')
		else:
			sys.stdout.write(odd_row_string + '\n')
			sys.stdout.write(even_row_string + '\n')
		sys.stdout.write('\n')
	sys.stdout.write('\n\n')

	return






######
# This function inputs the node array and the number of 42bp zones
# and returns a token pointer pair array
######
def token_pointer_pair_array(node_ra, num_zones, periodic_structure_flag):
	# Initialize the offset array
	even_offset_ra = [[ 0,  1], [ 0, -1], [-1, 0]]
	odd_offset_ra  = [[ 0, -1], [ 0,  1], [ 1, 0]]
	offset_ra = [even_offset_ra, odd_offset_ra]

			
	#Initialize the token pointer pair array
	strand_list_ra = []
	for sub_node_ra in node_ra:
		for node in sub_node_ra:
			if node != '.' and strand_list_ra.count(node) == 0:
				strand_list_ra.append(node)
	num_strands = len(strand_list_ra)

	TPP_ra = []
	num_subzones = num_zones*6
	for strand_num in range(num_strands):
		sub_TPP_ra = []
		for token_num in range(num_subzones):
			if strand_num%2 == 0:
				previous_TP = [strand_num, (token_num + 1)%num_subzones]
				if previous_TP[1] == 0 and periodic_structure_flag == False:
					previous_TP[1] = -1
				next_TP = [strand_num, (token_num - 1)%num_subzones]
				if next_TP[1] == (num_subzones - 1) and periodic_structure_flag == False:
					next_TP[1] = -1
				sub_TPP_ra.append([previous_TP, next_TP])
			else:
				previous_TP = [strand_num, (token_num - 1)%num_subzones]
				if previous_TP[1] == (num_subzones - 1) and periodic_structure_flag == False:
					previous_TP[1] = -1
				next_TP = [strand_num, (token_num + 1)%num_subzones]
				if next_TP[1] == 0 and periodic_structure_flag == False:
					next_TP[1] = -1
				sub_TPP_ra.append([previous_TP, next_TP])
		TPP_ra.append(sub_TPP_ra)


	# Introduce crossovers based on the node array
	for donor_y in range(len(node_ra)):
		for donor_x in range(len(node_ra[0])):
			donor_strand_num = node_ra[donor_y][donor_x]
			if donor_strand_num != '.':
				for position in range(3):
					parity = (donor_y + donor_x)%2
					acceptor_y = donor_y + offset_ra[parity][position][0]
					acceptor_x = donor_x + offset_ra[parity][position][1]
					acceptor_strand_num = node_ra[acceptor_y][acceptor_x]
					if acceptor_strand_num != '.':
						for zone_num in range(num_zones):
							subzone_num = zone_num*6 + position*2 + 1 - parity
							TPP_ra[donor_strand_num][subzone_num][1] = [acceptor_strand_num, subzone_num]
							TPP_ra[acceptor_strand_num][subzone_num][0] = [donor_strand_num, subzone_num]
							
	return TPP_ra






######
# This function inputs the token pointer pair array
# and returns an array of token_pointer_paths
######
def token_pointer_path_array(TPP_ra):
	sub_visits_ra = [0 for i in range(len(TPP_ra[0]))]
	visits_ra = [sub_visits_ra[:] for i in range(len(TPP_ra))]
	path_ra = []


	for strand_num in range(len(TPP_ra)):
		for subzone_num in range(len(TPP_ra[0])):
			if visits_ra[strand_num][subzone_num] == 0:	
				previous_TP = TPP_ra[strand_num][subzone_num][0]
				next_TP = TPP_ra[strand_num][subzone_num][1]
				num_visits = 100 * len(path_ra) + 1
				visits_ra[strand_num][subzone_num] = num_visits
				sub_path_ra = [[strand_num, subzone_num]]
							
				upstream_done = False
				while not upstream_done:
					if previous_TP[1] == -1:
						upstream_done = True
					elif visits_ra[previous_TP[0]][previous_TP[1]] > 0:
						upstream_done = True
					else:
						sub_path_ra.insert(0, previous_TP)
						num_visits += 1
						visits_ra[previous_TP[0]][previous_TP[1]] = num_visits
						previous_TP = TPP_ra[previous_TP[0]][previous_TP[1]][0]
			
				downstream_done = False
				while not downstream_done:
					if next_TP[1] == -1:
						downstream_done = True
					elif visits_ra[next_TP[0]][next_TP[1]] > 0:
						downstream_done = True
					else:
						sub_path_ra.append(next_TP)
						num_visits += 1
						visits_ra[next_TP[0]][next_TP[1]] = num_visits
						next_TP = TPP_ra[next_TP[0]][next_TP[1]][1]
			
				# Make sure that the path begins between subzones, not in the middle of a subzone
				if (sub_path_ra[0][0] + sub_path_ra[0][1])%2 == 0:
					sub_path_ra = sub_path_ra[1:] + sub_path_ra[:1]
				path_ra.append(sub_path_ra)



	num_tokens_visited = 0
	for sub_path_ra in path_ra:
		num_tokens_visited += len(sub_path_ra)

	sys.stdout.write('The number of tokens visited is ' + str(num_tokens_visited) + '.\n')
			

	path_length_ra = [0 for i in range(1100)]
	for sub_path_ra in path_ra:
		path_length_ra[len(sub_path_ra)] += 1

	for length in range(1100):
		if path_length_ra[length] != 0:
			sys.stdout.write('The number of paths with length ' + str(length) + ' is ' + str(path_length_ra[length]) + '.\n')
	
	return path_ra	






######
# This function inputs the token pointer path array
# and returns a list of oligos as lists of six token pointers
######
def oligo_token_pointer_array(path_ra):
	OTP_ra = []
	for sub_path_ra in path_ra:
		for oligo_num in range(len(sub_path_ra)/6):
			OTP_ra.append(sub_path_ra[oligo_num*6:oligo_num*6 + 6])
		# Take care of the remainders
		if len(sub_path_ra)%6 != 0:
			OTP_ra.append(sub_path_ra[-(len(sub_path_ra)%6):])
	return OTP_ra






######
# This function checks to make sure each token is represented once and only once
# in the oligo token pointer lists
######
def check_token_representation(OTP_ra, TPP_ra):
	CTP_ra = []
	for sub_OTP_ra in OTP_ra:
		for TP in sub_OTP_ra:
			CTP_ra.append(TP)

	num_strands = len(TPP_ra)
	num_subzones = len(TPP_ra[0])
	problems = False

	# Check each token for single presence
	for strand_num in range(num_strands):
		for token_num in range(num_subzones):
			if CTP_ra.count([strand_num, token_num]) == 0:
				sys.stdout.write(str([strand_num, token_num]) + ' not present.\n')
				problems = True
			elif CTP_ra.count([strand_num, token_num]) > 1:
				sys.stdout.write(str([strand_num, token_num]) + ' present more than once.\n')
				problems = True
	
	if problems == False:
		sys.stdout.write('Each token is represented once and only once in the oligo token pointer lists.\n')

	return
	





######
# This function prints the oligo path on the strand token lattice
# It is fun to see how the paths twist around the lattice
# You can use this function to help debug your program
######
def print_path(sub_path_ra, TPP_ra):
	sys.stdout.write('\nONE PATH ARRAY\n')

	num_strands = len(TPP_ra)
	num_subzones = len(TPP_ra[0])
	num_path_tokens = len(sub_path_ra)
	
	# Initialize strand token lattice
	sub_token_visit_ra = ['.' for subzone_num in range(num_subzones)]
	token_visit_ra = [sub_token_visit_ra[:] for strand_num in range(num_strands)]


	# Assign visits
	for path_token_num in range(num_path_tokens):
		token = sub_path_ra[path_token_num]
		strand = token[0]
		subzone = token[1]
		token_visit_ra[strand][subzone] = path_token_num

	# Print out strand token lattice
	spacer = '   '
	for strand_num in range(num_strands):
		for subzone_num in range(num_subzones):
			visitor_string = str(token_visit_ra[strand_num][subzone_num])
			sys.stdout.write(visitor_string)
			sys.stdout.write(spacer[:4 - len(visitor_string)])
		sys.stdout.write('\n')
	
	sys.stdout.write('\n')
	
	return
	
# The idea here is to have a function that adds the numbers of one oligo path 
# to the appropriate places in the big grid array. Eventually this will be printed
# in main. Also it needs to be initialized in main. Oligo_path is the path of 
# one oligo, while grid_ra is the grid that is constantly being updated until 
# it is printed in main. oligo_num is number that will be inputed to the grid_ra.

def generate_oligo_path(oligo_path, oligo_num, grid_ra):
	num_path_tokens = len(oligo_path)
        
# Assign visits
        for path_token_num in range(num_path_tokens):
                token = oligo_path[path_token_num]
                strand = token[0]
                subzone = token[1]
                grid_ra[strand][subzone] = oligo_num
                
        
	return grid_ra


def print_all_oligos(grid_ra, num_strands, num_subzones):
        spacer = '   '
        for strand_num in range(num_strands):
                for subzone_num in range(num_subzones):
                        visitor_string = str(grid_ra[strand_num][subzone_num])
                        sys.stdout.write(visitor_string)
                        sys.stdout.write(spacer[:4 - len(visitor_string)])
                sys.stdout.write('\n')
        


####
# given an oligo to split, split it and return the new list of oligos
####
def split_oligo(new_OTP_ra, oligo_num, num_toks):

	print new_OTP_ra[oligo_num]

	original_oligo = new_OTP_ra[oligo_num]

	oligo_1 = original_oligo[:num_toks]
	oligo_2 = original_oligo[num_toks:]
	print oligo_1
	print'\n'
	print oligo_2
	
	new_OTP_ra[oligo_num] = oligo_1
	new_OTP_ra.insert(oligo_num + 1, oligo_2)
	return new_OTP_ra


sys.stdout.write('Honeycomb pointers module installed.\n')
Personal tools