import numpy as np import pandas import matplotlib.pyplot as pl import math import matplotlib as mpl from matplotlib import cm # Define some initial stuff and import analysis functions: #AA_key_old=['A','G','L','M','F','W','K','Q','E','S','P','V','I','C','Y','H','R','N','D','T'] AA_key=['A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'] # So we've got 46 orthogonal (or at least not super correlated) # dimensions. Add them in to the matrix # From "Hot spot prediction in protein-protein interactions by an ensemble system" # Liu et. al. BMC Systems Biology newnew=pandas.read_csv('app_data/new_props') oldold=pandas.read_csv('app_data/old_props') # Again, ugly to hard code in the number of properties (62) but # For now no harm no foul properties=np.zeros((62,20)) for i in np.arange(len(AA_key)): properties[0:16,i]=oldold[AA_key[i]] properties[16:,i]=newnew[AA_key[i]] AA_num_key_new=properties[1] AA_num_key=np.arange(20)+1 def get_sequence_dimension(re_poly): num_loops,num_clones=np.shape(re_poly) for i in np.arange(num_loops): if i == 0: max_len=len(re_poly[i,0]) elif i <= 5: max_len=np.vstack((max_len,len(re_poly[i,0]))) for j in np.arange(num_clones): if i == 0: if len(re_poly[i,j]) > max_len: max_len = len(re_poly[i,j]) else: if len(re_poly[i,j]) > max_len[i]: max_len[i] = len(re_poly[i,j]) max_len=max_len+3 # Add 3 here so there's a more pronounced space between loops ## SO NOW MAX LEN SHOULD HAVE THE MAXIMUM LENGTH OF EACH CDR LOOP ## if num_loops == 1: sequence_dim = max_len else: sequence_dim = int(sum(max_len)) return(max_len,sequence_dim) # NOTE, manuscript_arrange=False MUST be selected to run MHC analysis # I used to re-arrange the CDR loops for a more position-accurate # representation. In the current analysis, this isn't necessary. def gen_tcr_matrix(pre_poly,key=AA_num_key_new,binary=False, pre_mono=[],giveSize=[],return_Size=False,manuscript_arrange=False): # Do this so that if giveSize == []: if binary: # NEW ADDITION TO CLEAN THINGS UP A BIT # max_len1 = get_sequence_dimension(pre_poly)[0] max_len2 = get_sequence_dimension(pre_mono)[0] max_lenp=np.zeros(len(max_len1)) for i in np.arange(len(max_len1)): max_lenp[i]=int(max(max_len1[i],max_len2[i])) if type(max_lenp) == int: sequence_dim = max_lenp else: sequence_dim = int(sum(max_lenp)) else: max_lenp,sequence_dim=get_sequence_dimension(pre_poly) else: max_lenp = giveSize if type(max_lenp) == int: sequence_dim = max_lenp else: sequence_dim = int(sum(max_lenp)) final_poly=[] # initialize a variable # RE-ORGANIZE EVERYTHING SO IT LOOKS NICE IN MATRIX # But presumably, if you give a size, it's the size you want... if manuscript_arrange and giveSize == []: max_lenp = [max_lenp[0],max_lenp[1],max_lenp[2],max_lenp[5],max_lenp[4],max_lenp[3]] max_len = max_lenp for re_poly in [pre_poly,pre_mono]: # RE-ORGANIZE EVERYTHING SO IT LOOKS NICE IN MATRIX if manuscript_arrange: re_poly = [re_poly[0],re_poly[1],re_poly[2],re_poly[5],re_poly[4],re_poly[3]] numLoop,numClone = np.shape(re_poly) poly_PCA=np.zeros([numClone,sequence_dim]) for i in range(numClone): # For all of our polyreactive sequences... loop=0 for k in range(numLoop): # Scroll through all of the loops within a clone leng=len(re_poly[k][i]) #should be this if re-sampling # the int clause below is how we center align if type(max_len) == int: count = int((max_len-leng)/2.0) else: count=int((max_len[k]-leng)/2.0)+int(sum(max_len[:k])) for m in re_poly[k][i]: # SO IS THIS ONE for bootstrapping for j in range(len(key)): if m==AA_key[j]: poly_PCA[i][count]=key[j] count=count+1 loop=loop+1 if binary: # Unfortunate naming here for the binary case # but it is what it is... if final_poly == []: final_poly = poly_PCA else: final_mono = poly_PCA else: break if binary and return_Size: return(final_poly,final_mono,max_lenp) elif binary: return(final_poly,final_mono) elif return_Size: return(poly_PCA,max_lenp) else: return(poly_PCA) def calculate_shannon(poly_PCA): clones,aas = np.shape(poly_PCA) prob_poly_full=np.zeros((clones,aas,21)) # We technically have 21 entries, 20 AAs (1-20) and spaces 0... Take a look at all of that AAs=np.arange(0,21) #print(AAs) for i in np.arange(clones): for j in np.arange((aas)): for k in AAs: if poly_PCA[i,j]==k: prob_poly_full[i,j,k]=prob_poly_full[i,j,k]+1 poly_count=np.sum(prob_poly_full,axis=0)/clones shannon_poly=np.zeros(len(poly_count)) for i in np.arange(len(poly_count)): for j in np.arange(len(poly_count[0])): if poly_count[i,j]==0: continue shannon_poly[i]=shannon_poly[i]+(-poly_count[i,j]*math.log(poly_count[i,j],2)) return(shannon_poly,poly_count) def calculate_MI(poly_PCA): shannon_poly,poly_count=calculate_shannon(poly_PCA) # Need to start removing explicit "126" sized matrices clones,aas = np.shape(poly_PCA) AAs=np.arange(0,21) MI_final_poly=np.zeros((aas,aas)) poly_count_cond=np.zeros((21,aas,aas,21)) save_count=np.zeros((21,aas)) for location in np.arange(aas): # Copy same stuff as above, but now also calculate a conditional (cond) probability cond_count_poly=0 #MI_final_poly=np.zeros((20,126)) conditional_entropy_poly=np.zeros((21,aas)) for res in AAs: prob_poly_cond=np.zeros((clones,aas,21)) cond_count_poly=0 for i in np.arange(clones): # Cycle through all clones if poly_PCA[i,location]==res: # Of those clones, only look at the ones that have AA 15 at position 73 cond_count_poly=cond_count_poly+1 for j in np.arange((aas)): # Then scroll through all of the positions to get conditional prob for k in AAs: if poly_PCA[i,j]==k: prob_poly_cond[i,j,k]=prob_poly_cond[i,j,k]+1 #So the matrix is built, now let's do our entropy... # for every position, we have a count... Can check and averaging over each position does add to 1. #if location == 75 and res == 5: # return(prob_poly_cond,cond_count_poly) if cond_count_poly == 0: continue poly_count_cond[res,location,:,:]=np.sum(prob_poly_cond,axis=0)/cond_count_poly save_count[res,location]=cond_count_poly MI_poly=np.zeros(len(poly_count_cond[res,location,:,:])) # Note, really should change this... MI poly is a misnomer here. This is one of the conditional entropy terms... for i in np.arange(len(poly_count_cond[res,location,:,:])): for j in np.arange(len(poly_count_cond[res,location,:,:][0])): if poly_count_cond[res,location,:,:][i,j]==0: continue MI_poly[i]=MI_poly[i]+(-poly_count_cond[res,location,:,:][i,j]*math.log(poly_count_cond[res,location,:,:][i,j],2)) condition_poly=np.zeros(len(shannon_poly)) condition_poly[location]=shannon_poly[location] #if cond_count_poly==0: # continue ############################## NOTE, I THINK THIS IS WHERE MY ISSUE IS ##################### # I need to sum up my MI_poly (which again, is actually conditional entropy) and THEN subtract from the standard entropy... conditional_entropy_poly[res]=MI_poly # Land at Shannon Entropy - global probability of a given residue * the conditional shannon entropy MI_final_poly[location]=shannon_poly-np.matmul(poly_count[location],conditional_entropy_poly)#-condition_poly return(MI_final_poly,poly_count_cond,poly_count) def joint_prob(poly_PCA): #Joint Probability Distribution Time! clones,aas = np.shape(poly_PCA) # Might be a bit of a beast to calculate... joint_probs = np.zeros((aas,21,aas,21)) # Should be position x number amino acids x AAs=np.arange(0,21) for cloneX in np.arange(clones): # Go through each clone in the matrix for posX in np.arange(aas): # GIVEN this first position in the matrix for resX in AAs: # cycle through all of the amino acids if poly_PCA[cloneX,posX] == resX: # test for existence to speed things up a bit for posY in np.arange(aas): # cycle through every position for resY in AAs: #cycle through every residue at these positions if poly_PCA[cloneX,posY] == resY: joint_probs[posX,resX,posY,resY] = joint_probs[posX,resX,posY,resY] + 1 joint=joint_probs/clones # divide out by the number of clones and number of Y positions return(joint) def gen_dset_props(poly_PCA,props=properties,stdev=False): AA_num_key=np.arange(20)+1 clones,aas = np.shape(poly_PCA) poly_prop=np.zeros([len(props),aas]) for i in np.arange(len(props)): # For all of our properties... for j in np.arange(clones): # for every clone for k in np.arange(aas): # for every position for m in AA_num_key: if poly_PCA[j,k]==m: poly_prop[i,k]=poly_prop[i,k]+props[i][m-1] poly_prop=poly_prop/clones if stdev: poly_prop_stdev=np.zeros([len(props),aas]) for i in np.arange(len(props)): # For all of our properties... for j in np.arange(clones): # for every clone for k in np.arange(aas): # for every position for m in AA_num_key: if poly_PCA[j,k]==m: #stdev = sqrt([sum(x-xavg)**2]/[N-1])... So below term is computing that sum poly_prop_stdev[i,k]=poly_prop_stdev[i,k]+(props[i][m-1]-poly_prop[i,k])**2 poly_prop_stdev=np.sqrt(poly_prop_stdev/(clones-1)) return(poly_prop,poly_prop_stdev) else: return(poly_prop) def gen_clone_props(poly_PCA): # Difference between this and the above function is how the # average is being taken... Either average over clones or average # over positions within the clones. clones,aas = np.shape(poly_PCA) # This whole code block used to be outside the function, but kept running into a very weird issue # where it seemed like the "props" were changing in a strange way. properties=np.zeros((62,20)) for i in np.arange(len(AA_key)): properties[0:16,i]=oldold[AA_key[i]] properties[16:,i]=newnew[AA_key[i]] props = properties[1:] # Re-normalize the properties for use in the matrix... for i in np.arange(len(props)): props[i] = props[i]/np.linalg.norm(props[i]) poly_prop_pca=np.zeros([len(props),clones]) for i in np.arange(len(props)): # For all of our properties... for j in np.arange(clones): # for every clone for k in np.arange(aas): # for every position for m in AA_num_key: if poly_PCA[j,k]==m: poly_prop_pca[i,j]=poly_prop_pca[i,j]+props[i,m-1] poly_prop_pca=poly_prop_pca/aas # 30 is a placeholder for now (30 positions tested) #single_pos_charge=np.zeros([30,clones]) #single_pos_phob=np.zeros([30,clones]) #for j in np.arange(clones): # for every clone # a=0 # a is here as a placeholder # for k in np.arange(68,98): # for every position # for m in AA_num_key: # if poly_PCA[j,k]==m: # So I think that 0 and 1 should correspond to charge and hydrophobicity... # single_pos_charge[a,j]=props[1,m-1] # single_pos_phob[a,j]=AA_num_key_new[m-1] # a=a+1 return(poly_prop_pca)#,single_pos_charge,single_pos_phob) def get_props(): # Identical to what we've got above, but it lets us also load it in to # the current variable space... AA_key=['A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V'] # So we've got 46 orthogonal (or at least not super correlated) # dimensions. Add them in to the matrix # From "Hot spot prediction in protein-protein interactions by an ensemble system" # Liu et. al. BMC Systems Biology (these are "new props") newnew=pandas.read_csv('app_data/new_props') oldold=pandas.read_csv('app_data/old_props') props=np.zeros((62,20)) for i in np.arange(len(AA_key)): props[0:16,i]=oldold[AA_key[i]] props[16:,i]=newnew[AA_key[i]] AA_num_key=np.arange(20)+1 AA_num_key_new=props[1] return(AA_key,AA_num_key,AA_num_key_new,props) def prop_patterning(mono_PCA,poly_PCA,mat_size=100,props=properties[1:],ridZero=False,win_size = 3,returnBig=False): # Try to maximize differences across the properties by looking at patterning... # Re-normalize the properties for use in the matrix... for i in np.arange(len(props)): props[i] = props[i]-np.average(props[i]) props[i] = props[i]/np.linalg.norm(props[i]) # Since we'll be averaging shit, let's get rid of all the zeros... # However, that's going to result in bleed-over across the loops... Is this good or bad? # This is also going to have a strange effect based upon loop length... # WHATEVER, Try both if ridZero: poly_pca_NEW=np.transpose(poly_PCA)[~np.all(np.transpose(poly_PCA) == 0,axis=1)] poly_pca_NEW=np.transpose(poly_pca_NEW) mono_pca_NEW=np.transpose(mono_PCA)[~np.all(np.transpose(mono_PCA) == 0,axis=1)] mono_pca_NEW=np.transpose(mono_pca_NEW) else: poly_pca_NEW = poly_PCA mono_pca_NEW = mono_PCA poly_dim1,poly_dim2=np.shape(poly_pca_NEW) mono_dim1,mono_dim2=np.shape(mono_pca_NEW) # So this is where we should be able to do the averaging poly_prop_masks=np.zeros([len(props),poly_dim1,int(poly_dim2/win_size)]) mono_prop_masks=np.zeros([len(props),mono_dim1,int(mono_dim2/win_size)]) for i in np.arange(len(props)): # For all of our properties... for j in np.arange(poly_dim1): # for every clone for k in np.arange(int(poly_dim2/win_size)): # for every position for win_slide in np.arange(win_size): # Make sure we don't go too far... if k*win_slide+win_slide > int(poly_dim2/win_size): continue for m in AA_num_key: if poly_pca_NEW[j,k*win_size+win_slide]==m: poly_prop_masks[i,j,k]=poly_prop_masks[i,j,k]+props[i,m-1] for i in np.arange(len(props)): # For all of our properties... for j in np.arange(mono_dim1): # for every clone for k in np.arange(int(mono_dim2/win_size)): # for every position for win_slide in np.arange(win_size): # Make sure we don't go too far... if k*win_slide+win_slide > int(mono_dim2/win_size): continue for m in AA_num_key: if mono_pca_NEW[j,k*win_size+win_slide]==m: mono_prop_masks[i,j,k]=mono_prop_masks[i,j,k]+props[i,m-1] # And now we actually have to do some sort of intelligent analysis on these masks... mono_prop_line = np.average(mono_prop_masks/win_size,axis=1) poly_prop_line = np.average(poly_prop_masks/win_size,axis=1) line_diff = poly_prop_line - mono_prop_line # Pull out the biggest differences in these properties... max_diffs = np.zeros((mat_size,3)) for i in np.arange(len(line_diff)): for j in np.arange(len(line_diff[0])): if line_diff[i,j] == 0: continue elif abs(line_diff[i,j]) > min(abs(max_diffs[:,0])): # Find the absolute value minimum. new_pos=np.where(abs(max_diffs[:,0]) == min(abs(max_diffs[:,0]))) # Put_it_here should get the location of the min... # It should always be in max_diffs[X,0] put_it_here=int(new_pos[0][0]) max_diffs[put_it_here,0] = line_diff[i,j] max_diffs[put_it_here,1] = i # Need the *win_size to revert back to physical location on the CDR loops... # Actually NOT true if we get rid of the zeros... little trickier there, come back to it max_diffs[put_it_here,2] = j*win_size # Now bring on back these properties for classification... new_mat_poly=np.zeros((mat_size,len(poly_pca_NEW))) new_mat_mono=np.zeros((mat_size,len(mono_pca_NEW))) for j in np.arange(len(poly_pca_NEW)): # for every clone for k in np.arange(len(max_diffs)): for win_slide in np.arange(win_size): # This really shouldn't happen, but just in case... if max_diffs[k,2]+win_slide > len(poly_pca_NEW[0]): continue for m in AA_num_key: if poly_pca_NEW[j,int(max_diffs[k,2]+win_slide)]==m: # So I think that 0 and 1 should correspond to charge and hydrophobicity... new_mat_poly[k,j]=new_mat_poly[k,j] + props[int(max_diffs[k,1]),m-1] for j in np.arange(len(mono_pca_NEW)): # for every clone for k in np.arange(len(max_diffs)): for win_slide in np.arange(win_size): # This really shouldn't happen, but just in case... if max_diffs[k,2]+win_slide > len(mono_pca_NEW[0]): continue for m in AA_num_key: if mono_pca_NEW[j,int(max_diffs[k,2]+win_slide)]==m: # So I think that 0 and 1 should correspond to charge and hydrophobicity... new_mat_mono[k,j]=new_mat_mono[k,j] + props[int(max_diffs[k,1]),m-1] if returnBig: return(new_mat_mono/win_size,new_mat_poly/win_size,max_diffs,poly_prop_masks,mono_prop_masks) else: return(new_mat_mono/win_size,new_mat_poly/win_size,max_diffs) # Alright so this below section is identical to the above code... For now. Need to rationalize how to pair/score def prop_pairing(ALL_mono,ALL_poly,mat_size=100,props=properties[1:],win_size = 3): # Try to maximize differences across the properties by looking at patterning... max_len1=get_sequence_dimension(ALL_poly)[0] max_len2=get_sequence_dimension(ALL_mono)[0] max_lenp=np.zeros(6) for i in np.arange(6): max_lenp[i]=max(max_len1[i],max_len2[i]) max_size=int(max(max_lenp)) # TEMPORARY #max_size=21 poly_PCA=gen_tcr_matrixOLD(ALL_poly,max_size,key=AA_num_key_new) mono_PCA=gen_tcr_matrixOLD(ALL_mono,max_size,key=AA_num_key_new) # For this one (and maybe the one above) we should reshape the matrix poly_dim1,poly_dim2=np.shape(poly_PCA) mono_dim1,mono_dim2=np.shape(mono_PCA) # This reshapes the matrices in to #clones by CDR loop by arbitrary position poly_pca_re = poly_PCA.reshape(int(poly_dim1),6,int(poly_dim2/6)) mono_pca_re = mono_PCA.reshape(int(mono_dim1),6,int(mono_dim2/6)) # FOR THIS SPECIFIC APPLICATION, I THINK ONLY PROPERTIES 1,2 apply (charge and hydrophobicity) props=props[1:3,:] # So this is where we should be able to do the averaging... So now win_size can't be over 21 poly_prop_masks=np.zeros([len(props),int(poly_dim1),6,int(poly_dim2/6/win_size)]) mono_prop_masks=np.zeros([len(props),int(mono_dim1),6,int(mono_dim2/6/win_size)]) # Make the property mask for the polyreactive matrix for i in np.arange(len(props)): # For all of our properties... for j in np.arange(poly_dim1): # for every clone for k in np.arange(int(poly_dim2/6/win_size)): # for every position for loop in np.arange(6): # for all 6 loops for win_slide in np.arange(win_size): # Make sure we don't go too far... if k*win_slide+win_slide > int(poly_dim2/6/win_size): continue for m in AA_num_key: if poly_pca_re[j,loop,k*win_size+win_slide]==m: poly_prop_masks[i,j,loop,k]=poly_prop_masks[i,j,loop,k]+props[i,m-1] # Make the property mask for the monoreactive matrix for i in np.arange(len(props)): # For all of our properties... for j in np.arange(mono_dim1): # for every clone for k in np.arange(int(mono_dim2/6/win_size)): # for every position for loop in np.arange(6): # for all 6 loops for win_slide in np.arange(win_size): # Make sure we don't go too far... if k*win_slide+win_slide > int(mono_dim2/6/win_size): continue for m in AA_num_key: if mono_pca_re[j,loop,k*win_size+win_slide]==m: mono_prop_masks[i,j,loop,k]=mono_prop_masks[i,j,loop,k]+props[i,m-1] # Alright now I've actually got to go through and determine if there are positive interactions between loops # FOR NOW, make the assumption that only AAs at the same position "see" each other... This is less egregious when averaging over windows Mmask_dim1,Mmask_dim2,Mmask_dim3,Mmask_dim4 = np.shape(mono_prop_masks) Pmask_dim1,Pmask_dim2,Pmask_dim3,Pmask_dim4 = np.shape(poly_prop_masks) poly_pair_score=np.zeros([6,6,Pmask_dim1,Pmask_dim2,Pmask_dim4]) mono_pair_score=np.zeros([6,6,Mmask_dim1,Mmask_dim2,Mmask_dim4]) for i in np.arange(Mmask_dim1): # For all of our properties for j in np.arange(Mmask_dim2): # for every clone for k in np.arange(Mmask_dim4): # for every position for loop1 in np.arange(6): # This and the bottom loop should be enough to find out loop-loop scoring for loop2 in np.arange(loop1+1): # The plus 1 here also includes self-loop if i == 0: # i=0 is charge as a property, so like-like is bad... if (mono_prop_masks[i,j,loop1,k] > 0.0) & (mono_prop_masks[i,j,loop2,k] > 0.0): mono_pair_score[loop1,loop2,i,j,k] = -1 mono_pair_score[loop2,loop1,i,j,k] = -1 elif (mono_prop_masks[i,j,loop1,k] < 0.0) & (mono_prop_masks[i,j,loop2,k] < 0.0): mono_pair_score[loop1,loop2,i,j,k] = -1 mono_pair_score[loop2,loop1,i,j,k] = -1 elif (mono_prop_masks[i,j,loop1,k] < 0.0) & (mono_prop_masks[i,j,loop2,k] > 0.0): mono_pair_score[loop1,loop2,i,j,k] = 1 mono_pair_score[loop2,loop1,i,j,k] = 1 elif (mono_prop_masks[i,j,loop1,k] > 0.0) & (mono_prop_masks[i,j,loop2,k] < 0.0): mono_pair_score[loop1,loop2,i,j,k] = 1 mono_pair_score[loop2,loop1,i,j,k] = 1 elif i == 1: # i=1 is hydrophobicity as a property, so like-like is good! if (mono_prop_masks[i,j,loop1,k] > 0.0) & (mono_prop_masks[i,j,loop2,k] > 0.0): mono_pair_score[loop1,loop2,i,j,k] = 1 mono_pair_score[loop2,loop1,i,j,k] = 1 elif (mono_prop_masks[i,j,loop1,k] < 0.0) & (mono_prop_masks[i,j,loop2,k] < 0.0): mono_pair_score[loop1,loop2,i,j,k] = 1 mono_pair_score[loop2,loop1,i,j,k] = 1 elif (mono_prop_masks[i,j,loop1,k] < 0.0) & (mono_prop_masks[i,j,loop2,k] > 0.0): mono_pair_score[loop1,loop2,i,j,k] = -1 mono_pair_score[loop2,loop1,i,j,k] = -1 elif (mono_prop_masks[i,j,loop1,k] > 0.0) & (mono_prop_masks[i,j,loop2,k] < 0.0): mono_pair_score[loop1,loop2,i,j,k] = -1 mono_pair_score[loop2,loop1,i,j,k] = -1 for i in np.arange(Pmask_dim1): # For all of our properties for j in np.arange(Pmask_dim2): # for every clone for k in np.arange(Pmask_dim4): # for every position for loop1 in np.arange(6): # This and the bottom loop should be enough to find out loop-loop scoring for loop2 in np.arange(loop1+1): # The plus 1 here also includes self-loop if i == 0: # i=0 is charge as a property, so like-like is bad... if (poly_prop_masks[i,j,loop1,k] > 0.0) & (poly_prop_masks[i,j,loop2,k] > 0.0): poly_pair_score[loop1,loop2,i,j,k] = -1 poly_pair_score[loop2,loop1,i,j,k] = -1 elif (poly_prop_masks[i,j,loop1,k] < 0.0) & (poly_prop_masks[i,j,loop2,k] < 0.0): poly_pair_score[loop1,loop2,i,j,k] = -1 poly_pair_score[loop2,loop1,i,j,k] = -1 elif (poly_prop_masks[i,j,loop1,k] < 0.0) & (poly_prop_masks[i,j,loop2,k] > 0.0): poly_pair_score[loop1,loop2,i,j,k] = 1 poly_pair_score[loop2,loop1,i,j,k] = 1 elif (poly_prop_masks[i,j,loop1,k] > 0.0) & (poly_prop_masks[i,j,loop2,k] < 0.0): poly_pair_score[loop1,loop2,i,j,k] = 1 poly_pair_score[loop2,loop1,i,j,k] = 1 elif i == 1: # i=1 is hydrophobicity as a property, so like-like is good! if (poly_prop_masks[i,j,loop1,k] > 0.0) & (poly_prop_masks[i,j,loop2,k] > 0.0): poly_pair_score[loop1,loop2,i,j,k] = 1 poly_pair_score[loop2,loop1,i,j,k] = 1 elif (poly_prop_masks[i,j,loop1,k] < 0.0) & (poly_prop_masks[i,j,loop2,k] < 0.0): poly_pair_score[loop1,loop2,i,j,k] = 1 poly_pair_score[loop2,loop1,i,j,k] = 1 elif (poly_prop_masks[i,j,loop1,k] < 0.0) & (poly_prop_masks[i,j,loop2,k] > 0.0): poly_pair_score[loop1,loop2,i,j,k] = -1 poly_pair_score[loop2,loop1,i,j,k] = -1 elif (poly_prop_masks[i,j,loop1,k] > 0.0) & (poly_prop_masks[i,j,loop2,k] < 0.0): poly_pair_score[loop1,loop2,i,j,k] = -1 poly_pair_score[loop2,loop1,i,j,k] = -1 # And now we actually have to do some sort of intelligent analysis on these masks... #mono_prop_line = np.average(mono_prop_masks/win_size,axis=1) #poly_prop_line = np.average(poly_prop_masks/win_size,axis=1) return(poly_pair_score,mono_pair_score,poly_prop_masks,mono_prop_masks) # OK, SO RATHER THAN change up this entire script, lemme just bring back in the OLD # gen_tcr_matrix script and use the old prop pairing... # What do I mean by this... It means one of the scripts requires each loop to be # the same number of entries. New way of making tcr_matrix doesn't do that. # So here, we just make 6 loop entries that are of length = max_len def gen_tcr_matrixOLD(re_poly,max_len,key=AA_num_key_new): poly_PCA=np.zeros([len(re_poly[0]),6*max_len]) for i in range(len(re_poly[:][0])): # For all of our polyreactive sequences... loop=0 for k in [0,1,2,5,4,3]: # Scroll through all of the loops within a clone #if poly[k][i][0]=='': # continue #leng=len(re_poly[k][i][0]) # THIS LINE IS AN ISSUE WHEN RESAMPLING leng=len(re_poly[k][i]) #should be this if re-sampling count=0+max_len*loop if leng min(abs(max_diffs[:,0])): # Find the absolute value minimum. new_pos=np.where(abs(max_diffs[:,0]) == min(abs(max_diffs[:,0]))) # Put_it_here should get the location of the min... # It should always be in max_diffs[X,0] put_it_here=int(new_pos[0][0]) max_diffs[put_it_here,0] = line_diff[i] max_diffs[put_it_here,1] = i # WE NO LONGER CARE ABOUT REVERTING TO PHYSICAL LOCATION (for this version) # Need the *win_size to revert back to physical location on the CDR loops... #max_diffs[put_it_here,2] = i*win_size #Think I can just delete all of the other shit and return max_diffs # Then I can just pull out those locations. return(max_diffs)