<jupyter_start><jupyter_text>STEP 1: Load Beacon, MAF, Reference and other cached variables<jupyter_code># CEU
beacon = pd.read_csv(join(ceuPath, "Beacon_164.txt"), index_col=0, delim_whitespace=True)
maf = pd.read_csv(join(ceuPath, "MAF.txt"), index_col=0, delim_whitespace=True)
reference = pickle.load(open(join(ceuPath, "reference.pickle"),"rb"))
binary = np.logical_and(beacon.values != reference, beacon.values != "NN").astype(int)
ternary = binary.copy()
ternary[beacon.values=="NN"] = -1

maf.rename(columns = {'referenceAllele':'major', 'referenceAlleleFrequency':'major_freq', 'otherAllele':'minor', 'otherAlleleFrequency':'minor_freq'}, inplace = True)

beaconPeople = np.arange(65)
otherPeople = np.arange(99)+65
allPeople = np.arange(164)
ternary = np.zeros(beacon.shape, dtype=int)
ternary[beacon != reference] = 1
ternary[beacon == "NN"] = -1
pca = PCA(n_components=4)
tr1 = pca.fit_transform(ternary.T)
plt.scatter(tr1[:, 0], tr1[:, 1], alpha=0.4)
kmeans = KMeans(3)
y_kmeans = kmeans.fit_predict(tr1)
plt.scatter(tr1[:, 0], tr1[:, 1], c=y_kmeans, s=50, cmap='viridis')
centers = kmeans.cluster_centers_
print(Counter(y_kmeans))
print(centers)

# Get indices of homogeneous groups
g1 = allPeople[np.where(y_kmeans == 0)[0]]
g2 = allPeople[np.where(y_kmeans == 1)[0]]
g3 = allPeople[np.where(y_kmeans == 2)[0]]
print(g1.shape, g2.shape, g3.shape)<jupyter_output>Counter({0: 81, 1: 75, 2: 8})
[[-976.14 11.53 -0.02 -0.02]
 [938.20 -171.92 0.20 0.06]
 [1087.75 1495.00 -1.68 -0.29]]
(81,) (75,) (8,)<jupyter_text>Chromosome Seperation<jupyter_code>chromosome_index = 1
ind = (maf["chromosome"] == "chr"+str(chromosome_index)).values
print(np.sum(ind)," SNP's exist in chromosome ", chromosome_index)

reference = reference[ind]
beacon = beacon.loc[ind]
extra = extra.loc[ind]
giant = giant.loc[ind]
maf = maf.loc[ind]<jupyter_output><empty_output><jupyter_text>STEP 1.2: Function Definitions<jupyter_code># Beacon operations
def queryBeacon(beacon_people):
    return binary[:, beacon_people].any(axis=1)

def getMutationAt(index):
    temp = maf.iloc[index]
    if temp["minor_freq"] == temp["maf"]:
        return temp["minor"] + temp["minor"] 
    else:
        return temp["major"] + temp["major"] 

def div(n, d):
    return n / d if d else 0

def rpaCalculate(tp,fp,tn,fn):
    recall = div(tp,(tp+fn)) 
    precision = div(tp,(tp+fp))
    accuracy = div((tp+tn),(tp+fp+tn+fn))
    return recall, precision, accuracy

# Performance method
def performance(person, reconstruction, eval_pos, reference):
    ind = np.logical_and(person[eval_pos] != np.squeeze(reference)[eval_pos], person[eval_pos] != "NN")
    tp = np.sum(reconstruction[eval_pos][ind] != np.squeeze(reference)[eval_pos][ind])
    fn = np.sum(ind) - tp
    fp = np.sum(reconstruction[eval_pos][~ind] != np.squeeze(reference)[eval_pos][~ind])
    tn = np.sum(~ind) - fp
    return tp, fp, tn, fn

def performance_f(test_people, reconstructed, add_count, cluster_count, eval_pos):
    total_values = np.zeros((4))
    best_matches = []
    # For all people in victim set
    for i in range(add_count):
        all_combinations = np.zeros((4, cluster_count))
        rpa = np.zeros((3, cluster_count))
        # For each cluster obtained
        for j in range(cluster_count):
            all_combinations[:, j] = performance(test_people[i], reconstructed[j], eval_pos, reference)
            rpa[:, j] = rpaCalculate(*all_combinations[:, j])
        ind = np.argmax(rpa[0,:]*rpa[1,:])       #Best-match index
        best_matches.append(ind)
        total_values += all_combinations[:, ind] #Add total tp-fp-tn-fn
    recall, precision, accuracy = rpaCalculate(*total_values)
    print("Recall_Micro_Avg    =", round(recall, 2),"\nPrecision_Micro_Avg =", round(precision, 2))
    return (precision,recall,accuracy), total_values, best_matches<jupyter_output><empty_output><jupyter_text>STEP 2: Choose random people and send query to Beacon to obtain No-Yes answers<jupyter_code>def getNoYes(add_count, beacon_size):

    # Take people for added group
    added_people = otherPeople.copy()
    random.shuffle(added_people)
    added_people = added_people[:add_count]
    
    # Take people for beacon
    beacon_people = np.setdiff1d(allPeople, added_people)
    random.shuffle(beacon_people)
    beacon_people = beacon_people[:beacon_size]

    # Query Beacon initially
    before = queryBeacon(beacon_people)
    # Add people
    updated_beacon = np.concatenate([added_people,beacon_people])
    # Query Beacon again
    after = queryBeacon(updated_beacon)
    # Find No-Yes SNPs' indices
    no_yes_indices = np.where(np.logical_and(before==False, after==True))[0]
    yes_yes_indices = np.where(np.logical_and(before==True, after==True))[0]
    print("Number of No-Yes SNP's : ", len(no_yes_indices))
    
    return yes_yes_indices, no_yes_indices, added_people, beacon_people
def getNoYes2(add_count, beacon_size):

    # Take people for added group
    random.shuffle(g1)
    random.shuffle(g2)
    if add_count == 5:
        added_people = np.concatenate([g1[:3], g2[:2]])
    else:
        added_people = np.concatenate([g1[:(add_count // 2)], g2[:(add_count // 2)]])

    # Take people for beacon
    g1_ = np.setdiff1d(g1, added_people)
    random.shuffle(g1_)
    g2_ = np.setdiff1d(g2, added_people)
    random.shuffle(g2_)
    
    
    curBeacon = np.concatenate([g1_[:(beacon_size // 2)], g2_[:(beacon_size // 2)]])

    # Query Beacon initially
    before = queryBeacon(curBeacon)
    # Add people
    updatedBeacon = np.concatenate([added_people, curBeacon])
    # Query Beacon again
    after = queryBeacon(updatedBeacon)
    # Find No-Yes SNPs' indices
    no_yes_indices = np.where(np.logical_and(before == False, after == True))[0]
    yes_yes_indices = np.where(np.logical_and(before == True, after == True))[0]
    print("Number of No-Yes SNP's : ", len(no_yes_indices))

    return yes_yes_indices, no_yes_indices, added_people, curBeacon<jupyter_output><empty_output><jupyter_text>STEP 3: Correlation Model<jupyter_code>def builtSNPNetwork(no_yes_indices, model_ind, reference):
    model = ternary[no_yes_ind][:, model_ind].astype(float)
    model[model==-1] = np.nan
    x = pairwise_distances(model, metric = "sokalmichener", n_jobs=-1)
    x = 1-np.nan_to_num(x)
    return x
def baseline_method(no_yes_indices, add_count, cluster_count=None):
    c = maf.iloc[no_yes_indices]

    # Calculate probabilities of SNP possibilities
    greater = c.loc[c['major_freq'] >= c['minor_freq']]
    smaller = c.loc[c['major_freq'] < c['minor_freq']]

    greater["maj-maj"] = greater['major'] + "" + greater['major']
    greater["mean"] = pd.concat([greater['major'] + "" + greater['minor'], greater['minor'] + "" + greater['major']], axis=1).min(axis=1)
    greater["min-min"] = greater['minor'] + "" + greater['minor']
    greater["p1"] = greater['major_freq']**2
    greater["p2"] = 2*greater['major_freq']*greater['minor_freq']
    greater["p3"] = greater['minor_freq']**2

    smaller["maj-maj"] = smaller['minor'] + "" + smaller['minor']
    smaller["mean"] = pd.concat([smaller['major'] + "" + smaller['minor'], smaller['minor'] + "" + smaller['major']], axis=1).min(axis=1)
    smaller["min-min"] = smaller['major'] + "" + smaller['major']
    smaller["p1"] = smaller['minor_freq']**2
    smaller["p2"] = 2*smaller['major_freq']*smaller['minor_freq']
    smaller["p3"] = smaller['major_freq']**2

    tt = pd.concat([greater,smaller], axis=0)
    tt.sort_index(inplace=True)

    genome_possibilities = tt[["maj-maj", "mean", "min-min"]].values
    probabilities = tt[["p1","p2","p3"]].values

    mutations = tt[["mean", "min-min"]].values
    mutation_probs = tt[["p2","p3"]].values

    # Randomly reconstruct the people's genome
    bins = []
    cumulative = probabilities.cumsum(axis=1)
    for i in range(add_count):
        uniform = np.random.rand(len(cumulative), 1)
        choices = (uniform < cumulative).argmax(axis=1)
        reconstructed = np.choose(choices, genome_possibilities.T)
        bins.append(reconstructed)
    bins = np.array(bins)
    
    # Be sure that at least one person has the mutation
    equality = np.sum((bins == reference[no_yes_indices].T), axis=0)
    changed_indices = np.where(equality==add_count)[0]

    index_choices = np.random.randint(add_count, size=len(equality))[changed_indices]

    non_zeros = mutation_probs[np.sum(mutation_probs, axis=1) != 0]
    probs = (non_zeros.T / np.sum(non_zeros, axis=1).T).T

    zeros = np.zeros((mutation_probs.shape[0], 2))
    zeros[np.sum(mutation_probs, axis=1) != 0] = probs
    probs = zeros[changed_indices]

    cum = probs.cumsum(axis=1)
    uni = np.random.rand(len(cum), 1)
    choi = (uni < cum).argmax(axis=1)
    res = np.choose(choi, mutations[changed_indices].T)

    bins.T[changed_indices, index_choices] = res
    # Reconstruct
    reconstructed = np.array([reference.T[0] for i in range(add_count)])
    reconstructed.T[no_yes_indices] = bins.T
    return reconstructed<jupyter_output><empty_output><jupyter_text>Spectral Clustering<jupyter_code>def spectralClustering(no_yes_indices, add_count, x, reference, cluster_count=None):
    if not cluster_count:
        cluster_count = add_count
    sc = SpectralClustering(cluster_count, affinity='precomputed', n_init=100, n_jobs=-1)
    sc.fit(np.array(x))
    bins = []
    for i in range(cluster_count):
        temp = []
        for element in np.where(sc.labels_==i)[0]:
            temp.append(no_yes_indices[element])
        #print("Bin " + str(i) + " has " + str(len(temp)) + " SNP's")
        bins.append(temp)
    reconstructed = np.array([reference.T[0] for i in range(cluster_count)])
    for i in range(cluster_count):
        for j in bins[i]:
            reconstructed[i][j] = getMutationAt(j)
    return reconstructed<jupyter_output><empty_output><jupyter_text>Fuzzy Clustering<jupyter_code>def fuzzyClustering(no_yes_indices, add_count, x, reference, cluster_count=None):
    if not cluster_count:
        cluster_count = add_count
    fcm = FCM(n_clusters=cluster_count)
    fcm.fit(correlations)
    soft_clusters = fcm.u
    bins = [[] for i in range(cluster_count)]
    for i in range(len(soft_clusters)):
        maxPos = np.argmax(soft_clusters[i])
        if soft_clusters[i][maxPos] <= 0.5:
            for j in np.where(soft_clusters[i] > (soft_clusters[i][maxPos] * 2 / 3))[0]:
                bins[j].append(no_yes_indices[i])
        else:
            bins[maxPos].append(no_yes_indices[i])
    reconstructed = np.array([reference.T[0] for i in range(cluster_count)])
    for i in range(cluster_count):
        for j in bins[i]:
            reconstructed[i][j] = getMutationAt(j)
    return reconstructed<jupyter_output><empty_output><jupyter_text>E X P E R I M E N T<jupyter_code>add_count = 5
cluster_count = 5
beacon_size = 50
yes_yes_ind, no_yes_ind, added_people, beacon_people     = getNoYes2(add_count, beacon_size)<jupyter_output>Number of No-Yes SNP's :  14783<jupyter_text>All People<jupyter_code>model_ind = allPeople
correlations                              = builtSNPNetwork(no_yes_ind, model_ind, reference)
reconstructed_spectral                    = spectralClustering(no_yes_ind, add_count, correlations, reference)
(precision,recall,accuracy), _, matches   = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,cluster_count,no_yes_ind)<jupyter_output>Recall_Micro_Avg    = 0.75 
Precision_Micro_Avg = 0.72<jupyter_text>Without Added People<jupyter_code>model_ind = np.setdiff1d(allPeople, added_people)
model_ind = np.setdiff1d(model_ind, beacon_people)

correlations                              = builtSNPNetwork(no_yes_ind, model_ind, reference)
reconstructed_spectral                    = spectralClustering(no_yes_ind, add_count, correlations, reference)
(precision,recall,accuracy), _, matches   = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,cluster_count,no_yes_ind)<jupyter_output>Recall_Micro_Avg    = 0.55 
Precision_Micro_Avg = 0.42<jupyter_text>Genome Reconstruction<jupyter_code>add_count = 2
cluster_count = 2
beacon_size = 50
yes_yes_ind, no_yes_ind, added_people     = getNoYes(add_count, beacon_size)
model_ind = np.setdiff1d(otherPeople, added_people)
correlations                              = builtSNPNetwork(no_yes_ind, model_ind, reference)
reconstructed_spectral                    = spectralClustering(no_yes_ind, add_count, correlations, reference)
(precision,recall,accuracy), _, matches   = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,cluster_count,no_yes_ind)<jupyter_output><empty_output><jupyter_text>All Tests and Plots<jupyter_code>experiments = [(2,20),(3,30),(4,40),(5,50),(10,100),(20,100)]
res = []
for e in experiments:
    add_count = e[0]
    beacon_size = e[1]

    for i in range(20):
        yes_yes_ind, no_yes_ind, added_people = getNoYes(add_count, beacon_size)
        model_ind = np.setdiff1d(otherPeople, added_people)

        # Genome Reconstruction    
        correlations                              = builtSNPNetwork(no_yes_ind, model_ind, reference)
        reconstructed_spectral                    = spectralClustering(no_yes_ind, add_count, correlations, reference)
        (precision,recall,accuracy), _, matches   = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,add_count,no_yes_ind)<jupyter_output><empty_output><jupyter_text>Experiment 1: Vary Added People<jupyter_code>counts = [2,3,5,10,20]
beacon_size = 50
run_count = 20

results = np.zeros((3, len(counts), run_count, 3))

for i in range(len(counts)):
    for j in range(run_count):
        yes_yes_ind, no_yes_ind, added_people = getNoYes(counts[i], beacon_size)
        model_ind = np.setdiff1d(otherPeople, added_people)

        # Genome Reconstruction    
        correlations                = builtSNPNetwork(no_yes_ind, model_ind, reference)
        reconstructed_spectral      = spectralClustering(no_yes_ind, counts[i], correlations, reference)
        results[0, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,counts[i],counts[i],no_yes_ind) 

        # Baseline
        reconstructed_baseline      = baseline_method(no_yes_ind, counts[i])
        results[1, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_baseline,counts[i],counts[i],no_yes_ind)  

        # Fuzzy
        reconstructed_fuzzy         = fuzzyClustering(no_yes_ind, counts[i], correlations, reference)
        results[2, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_fuzzy, counts[i], counts[i], no_yes_ind)  
        
with open("../results/1C-Triple-VaryAdded.pickle", 'wb') as file:
    pickle.dump(results, file)<jupyter_output><empty_output><jupyter_text>Experiment 2: Vary Beacon Size<jupyter_code>counts = [25,50,75,100]
add_count = 5
run_count = 20

results = np.zeros((3, len(counts), run_count, 3))

for i in range(len(counts)):
    for j in range(run_count):
        yes_yes_ind, no_yes_ind, added_people = getNoYes(add_count, counts[i])
        model_ind = np.setdiff1d(otherPeople, added_people)

        # Genome Reconstruction    
        correlations                = builtSNPNetwork(no_yes_ind, model_ind, reference)
        reconstructed_spectral      = spectralClustering(no_yes_ind, add_count, correlations, reference)
        results[0, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,add_count,no_yes_ind) 

        # Baseline
        reconstructed_baseline      = baseline_method(no_yes_ind, add_count)
        results[1, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_baseline,add_count,add_count,no_yes_ind)  

        # Fuzzy
        reconstructed_fuzzy         = fuzzyClustering(no_yes_ind,add_count,correlations,reference)
        results[2, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_fuzzy,add_count,add_count,no_yes_ind)  
        
with open("../results/1C-Triple-VaryBeacon.pickle", 'wb') as file:
    pickle.dump(results, file)<jupyter_output><empty_output><jupyter_text>Experiment 3: Cluster Count<jupyter_code>counts = [1,2,3,4,5,10]
add_count = 5
run_count = 20
beacon_size = 50

results = np.zeros((3, len(counts), run_count, 3))

for i in range(len(counts)):
    for j in range(run_count):
        yes_yes_ind, no_yes_ind, added_people = getNoYes(add_count, beacon_size)
        model_ind = np.setdiff1d(otherPeople, added_people)

        # Genome Reconstruction    
        correlations                = builtSNPNetwork(no_yes_ind, model_ind, reference)
        reconstructed_spectral      = spectralClustering(no_yes_ind, add_count, correlations, reference)
        results[0, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T,reconstructed_spectral,add_count,counts[i],no_yes_ind) 

        # Baseline
        reconstructed_baseline      = baseline_method(no_yes_ind, add_count)
        results[1, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_baseline,add_count,counts[i],no_yes_ind)  

        # Fuzzy
        reconstructed_fuzzy         = fuzzyClustering(no_yes_ind,add_count,correlations,reference)
        results[2, i, j, :],_,_     = performance_f(beacon.iloc[:, added_people].values.T, reconstructed_fuzzy,add_count,counts[i],no_yes_ind)  
        
with open("../results/1C-Triple-VaryCluster.pickle", 'wb') as file:
    pickle.dump(results, file)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13649784238019533, 'fraction_numerical': 0.018112650465591643, 'mean_word_length': 4.2450863609291245, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '127291', 'license': 'MIT', 'n_tokens_mistral': 6625, 'n_tokens_neox': 6221, 'n_words': 1064, 'path': 'src/reconstruction_CEU.ipynb', 'repo_name': 'ciceklab/GenomeReconstruction', 'size': '17612'}