<jupyter_start><jupyter_text>IntroductionIn this example, the following stages are tested: 1- Loading File_1.zip and encoding it using a Reed-Solomon encoder. 2- Generating noisy reads from the encoded file. 3- Clustering the reads using the Locality-Sensitive Hashing (LSH) method. 4- Generating candidates from the clusters by performing multiple sequence alignment and majority voting. 5- Putting the candidates into a Reed-Solomon decoder and recovering the oiginal data.<jupyter_code>import numpy as np
from random import shuffle
from Bio.Align.Applications import MuscleCommandline
from Bio import AlignIO
from skbio.alignment import local_pairwise_align_ssw
from skbio import DNA
import itertools
import operator
import time
import zipfile
import subprocess as sp
import hashlib
from random import random<jupyter_output><empty_output><jupyter_text>EncodingNote that n, k, N, K, nus, and numblocks (this one is printed out when the following function is called) will be needed for decoding.<jupyter_code>def encode(infile):
    stat = sp.Popen(["./simulate/texttodna", "--n=16383", "--k=10977", "--N=20", "--K=18", "--nuss=6", "--l=4", \
                     "--encode", infile,"--output=./data/data_encoded.txt"],stdout = sp.PIPE)
    for line in stat.stdout.read().splitlines():
        print(line)
    return True
s = time.time()
encode("--input=./data/File_1.zip")
print("Encoding time:",time.time()-s)<jupyter_output>b'fieldpower: lookup tables generated.. '
b'fieldpower: lookup tables generated.. '
b'--------------------------------'
b'redundancy outer code: 0.492484 (= (n-k)/k)'
b'redundancy inner code: 0.111111 (= (N-K)/K)'
b'--------------------------------'
b'start encoding..'
b'infile:  ./data/File_1.zip'
b'outfile: ./data/data_encoded.txt'
b'Filesize in Byte: 99103'
b'encode block 0'
b'\t encode part  0'
b'\t encode part  1'
b'\t encode part  2'
b'\t encode part  3'
b'\t encode part  4'
b'\t encode part  5'
b'encoded block 0 of 1'
b'encoded 99103 Bytes to 1 blocks, resulting in 16383 DNA segments of length 20 each.'
Encoding time: 16.468868017196655<jupyter_text>Induce deletions and generate noisy reads (perturbing)<jupyter_code>def generate_noisy_reads(infile,p,n):
    # infile: encoded file
    # p: deletion probabilty per character
    # n: number of noisy reads per original sequence
    
    orig_seqs = []
    f = open(infile,"r")
    for seq in f.read().splitlines():
        orig_seqs += [seq]
    
    
    rand_replace = lambda c: c if random()<(1-p) else '' # this function induces deletions
    population = []
    for seq in orig_seqs:
        for i in range(n):
            s = ''.join([rand_replace(c) for c in seq])
            population += [s]
    shuffle(population)
    return orig_seqs,population
orig_seqs,reads = generate_noisy_reads("./data/data_encoded.txt",0.02,7)<jupyter_output><empty_output><jupyter_text>LSH clustering<jupyter_code>#===== assign numbers to shingles of each sequence=====#
def kmerDNA(seq,k=3):
    kmer = []
    for ell in range(len(seq)-k+1):
        nstr = seq[ell:ell+k]
        index = 0
        for j,c in enumerate(nstr):
            if c == 'A':
                i = 0
            elif c == 'C':
                i = 1
            elif c == 'G':
                i = 2
            elif c == 'T':
                i = 3
            else:
                index = -1
                break
            index += i*(4**j)
        kmer += [index]
    return kmer
#=====min-hash object=====#
class minhashsig():
    # min-hash of k-mers
    def __init__(self,m,k):
        # m is the number of signatures
        self.tables = [np.random.permutation(4**k) for i in range(m)]
        self.k = k
    def generate_signature(self,seq):
        kmer = kmerDNA(seq,self.k)
        sig = [ min([table[i] for i in kmer]) for table in self.tables]
        return sig
#=====pair detection=====#
def extract_similar_pairs(sigs,m,k_lsh,ell_lsh,maxsig):
    # sigs: minhash signatures
    # ell_lsh: number of LSH signatures
    # k_lsh: number of MH signatures to be concatenated
    # we use generatrs to yield a number of pairs at a time for the sake of memory efficiency
    
    pairs = set([])
    
    # generate ell_lsh random indices
    for ell in range(ell_lsh):
        pair_count = 0
        s = time.time()
        lshinds = np.random.permutation(m)[:k_lsh]
        # generate LSh signatures
        lshsigs = []
        for sig in sigs:
            lshsig = 0
            for i,lshind in enumerate(lshinds):
                lshsig += sig[lshind]*(maxsig**i)
            lshsigs += [lshsig]
        d = {}
        for ind,sig in enumerate(lshsigs):
            if sig in d:
                d[sig] += [ind]
            else:
                d[sig] = [ind]
        for candidates in d.values():
            cent = set([])
            if len(candidates) > 1:
                for pair in itertools.combinations(candidates,2):
                    cent.add(pair[0])
                    if len(cent)==1:
                        pairs.add(pair)
                    else:
                        break
                        
        yield pairs,ell
        pair_count += len(pairs)
        pairs = set([])
#=====form clusters based on pairs=====#
def center_cluster(pairs):
    clusters = {}
    hold = 0
    t_counter = 0
    ell_copy = 0
    pairsize = 0
    while not hold:
        
        try:
            out = next(pairs)
            pairs_sort = list(out[0])
            ell = out[1]
            pairsize += len(pairs_sort)
            pairs_sort.sort()
            s = time.time()
            for (u,v) in pairs_sort:
                if u in clusters:
                    clusters[u] += [v]
        
                if v in clusters:
                    clusters[v] += [u]
        
                if v not in clusters and u not in clusters:
                    clusters[u] = [v]
    
        except StopIteration:
            hold = 1
            print("clustering completed","---",pairsize,"pairs clustered")
        if ell==ell_copy:
            t_counter += time.time()-s
        else:
            print("Clustering time for LSH",ell_copy,":",t_counter,'\n')
            t_counter = time.time()-s
            ell_copy = ell
 
    return clusters

#=====LSH clustering (main function)=====#
def lsh_cluster(seqs,m,k,k_lsh=2,ell_lsh=4):
    # This is the main function
    maxsig = 4**k
    minhash = minhashsig(m,k)
    sigs = [minhash.generate_signature(seq[:40]) for seq in seqs]
    pairs = extract_similar_pairs(sigs,m,k_lsh,ell_lsh,maxsig)
    clusters = center_cluster(pairs)
    return clusters
# set up the parameters and call the lsh_cluster function
k_lsh = 4
sim = 0.5
ell_lsh = int(1/(sim**k_lsh))
m,k=50,5
start = time.time()
clusters = lsh_cluster(reads,m,k,k_lsh,ell_lsh)
end = time.time()

print("Runtime:",end-start)
len(clusters)<jupyter_output>Clustering time for LSH 0 : 0.0239560604095459 

Clustering time for LSH 1 : 0.029717683792114258 

Clustering time for LSH 2 : 0.030060768127441406 

Clustering time for LSH 3 : 0.04369521141052246 

Clustering time for LSH 4 : 0.03241896629333496 

Clustering time for LSH 5 : 0.03330802917480469 

Clustering time for LSH 6 : 0.03408098220825195 

Clustering time for LSH 7 : 0.03426837921142578 

Clustering time for LSH 8 : 0.03566384315490723 

Clustering time for LSH 9 : 0.0362248420715332 

Clustering time for LSH 10 : 0.03715825080871582 

Clustering time for LSH 11 : 0.035012245178222656 

Clustering time for LSH 12 : 0.0354459285736084 

Clustering time for LSH 13 : 0.03604912757873535 

Clustering time for LSH 14 : 0.03756213188171387 

clustering completed --- 1002366 pairs clustered
Runtime: 81.8294587135315<jupyter_text>Filtering the clustersIn this stage, we use the max_match function to filter the clusters. This function checks the simlarity of each cluster member and its cluster center based on local alignment.<jupyter_code># adding the center of each cluter to the cluster (also removing duplicates from each cluster)
clusts = [ [c] + list(set(clusters[c])) for c in clusters if len(clusters[c]) > 3 ]
#=====max matching=====# 
def max_match(seq1,seq2):
    # This function checks whether seq1 and seq2 are similar or not
    # Checking all pairs within a cluster dramatically increases the time complexity, 
    # so by default, in the next cell, we call this function to only check the pairs
    # that one of their members is the cluster center
    
    alignment,score,start_end_positions \
        = local_pairwise_align_ssw(DNA(seq1) , DNA(seq2) , match_score=2,mismatch_score=-3)
    a = str(alignment[0])
    b = str(alignment[1])
    ctr = 0
    for i,j in zip(a,b):
        if i==j:
            ctr += 1
    return ctr
th = 35 # filtering threshold

k = len(clusts)
s = time.time()
fclusts = []
for i,c in enumerate(clusts):
    cent = reads[c[0]]
    cc = [c[0]]
    for e in c[1:]:
        score = max_match(cent,reads[e])
        if score >= th:
            cc += [e]
    fclusts += [cc]
    if i%100 == 0:
        print(i/len(clusts))
print("filtering time for",k,"clusters:",time.time()-s)<jupyter_output>0.0
0.0026085141903171953
0.005217028380634391
0.007825542570951586
0.010434056761268781
0.013042570951585977
0.015651085141903172
0.018259599332220367
0.020868113522537562
0.023476627712854758
0.026085141903171953
0.02869365609348915
0.031302170283806344
0.03391068447412354
0.036519198664440734
0.03912771285475793
0.041736227045075125
0.044344741235392324
0.046953255425709516
0.049561769616026714
0.052170283806343906
0.054778797996661105
0.0573873121869783
0.059995826377295496
0.06260434056761269
0.06521285475792989
0.06782136894824708
0.07042988313856427
0.07303839732888147
0.07564691151919867
0.07825542570951587
0.08086393989983305
0.08347245409015025
0.08608096828046745
0.08868948247078465
0.09129799666110183
0.09390651085141903
0.09651502504173623
0.09912353923205343
0.10173205342237061
0.10434056761268781
0.10694908180300501
0.10955759599332221
0.1121661101836394
0.1147746243739566
0.11738313856427379
0.11999165275459099
0.12260016694490818
0.12520868113522537
0.12781719532554256[...]<jupyter_text>Aligning the clustersAs for the alignment stage, there are a number of options differing in accuracy and speed. For the sake of having the highest accuracy, we chose the Muscle function which is the most accurate multiple-alignment function among all python built-in functions.<jupyter_code>def multiple_alignment_muscle(cluster,out=False):
    # write cluster to file
    file = open("clm.fasta","w") 
    for i,c in enumerate(cluster):
        file.write(">S%d\n" % i)
        file.write(c)
        file.write("\n")
    file.close()

    muscle_exe = r"~/muscle3.8.31_i86linux64" # assuming you've already put this in the main directory
    output_alignment = "clmout.fasta"
    muscle_cline = MuscleCommandline(muscle_exe, input="clm.fasta", out=output_alignment)
    stdout, stderr = muscle_cline()
    msa = AlignIO.read(output_alignment, "fasta")
    if out:
        print(msa)
    alignedcluster = []
    for i in msa:
        alignedcluster += [i.seq]
    return alignedcluster
fresults = []
def align_clusters(clusters,masize = 15):
    ### align clusters, generate candidates
    for i, clusterinds in enumerate(clusters):
        cluster = [reads[i] for i in clusterinds]
        if len(cluster) < 2:
            continue
        
        ma = multiple_alignment_muscle(cluster[:masize])
        fresults.append(ma)
            
        if i % 100 == 0:
            print(i/len(clusters))
s = time.time()
align_clusters(fclusts,5)
print("Alignment Runtime:",time.time()-s)<jupyter_output>0.0
0.0026085141903171953
0.005217028380634391
0.007825542570951586
0.010434056761268781
0.013042570951585977
0.015651085141903172
0.018259599332220367
0.020868113522537562
0.023476627712854758
0.026085141903171953
0.02869365609348915
0.031302170283806344
0.03391068447412354
0.036519198664440734
0.03912771285475793
0.041736227045075125
0.044344741235392324
0.046953255425709516
0.049561769616026714
0.052170283806343906
0.054778797996661105
0.0573873121869783
0.059995826377295496
0.06260434056761269
0.06521285475792989
0.06782136894824708
0.07042988313856427
0.07303839732888147
0.07564691151919867
0.07825542570951587
0.08086393989983305
0.08347245409015025
0.08608096828046745
0.08868948247078465
0.09129799666110183
0.09390651085141903
0.09651502504173623
0.09912353923205343
0.10173205342237061
0.10434056761268781
0.10694908180300501
0.10955759599332221
0.1121661101836394
0.1147746243739566
0.11738313856427379
0.11999165275459099
0.12260016694490818
0.12520868113522537
0.12781719532554256[...]<jupyter_text>Majority mergingAt this stage, we merge each aligned cluster by putting up a voting for each position within the aligned sequences. Then, the fraction of all original sequences recovered is computed.<jupyter_code># This function returns the fraction of origignal squences recovered given a number of candidates
def fraction_recovered(candidates,orig_seqs):
    d = {}
    for seq in orig_seqs:
        d[seq] = 0
    for cand in candidates:
        if cand in d:
            d[cand] += 1
    av = sum([ d[seq]>0 for seq in d]) / len(d)
    print("Fraction of recovered sequences: ", av )
    if av>0:
        print("Fraction of recovered sequences: ", sum([ d[seq] for seq in d]) / len(d) / av )
def majority_merge(reads,weight = 0.4):
    # assume reads have the same length
    res = ""
    for i in range(len(reads[0])):
        counts = {'A':0,'C':0,'G':0,'T':0,'-':0,'N':0}
        for j in range(len(reads)):
            counts[reads[j][i]] +=1
        counts['-'] *= weight
        mv = max(counts.items(), key=operator.itemgetter(1))[0]
        if mv != '-':
            res += mv
    return res
candidates = []
for ma in fresults:
    candidates.append(majority_merge(ma,weight=0.5))
fraction_recovered( [seq[:60] for seq in candidates] , orig_seqs)
f = open("./data/data_drawnseg.txt","w")
for i,seq in enumerate(candidates):
    if len(seq) >= 60:
        f.write(seq[:60])
        f.write('\n')
f.close()<jupyter_output><empty_output><jupyter_text>Decoding<jupyter_code>def decode(infile):
    stat = sp.Popen(["./simulate/texttodna", "--decode",  "--n=16383", "--k=10977", "--N=20", "--K=18", "--nuss=6", "--l=4", \
                    "--numblocks=1", infile, "--output=./data/data_rec.zip"],stdout=sp.PIPE)
    for line in stat.stdout.read().splitlines():
        print(line)
    return True
s = time.time()
decode("--input=./data/data_drawnseg.txt")
print("Decoding time:",time.time()-s)<jupyter_output>b'fieldpower: lookup tables generated.. '
b'fieldpower: lookup tables generated.. '
b'--------------------------------'
b'redundancy outer code: 0.492484 (= (n-k)/k)'
b'redundancy inner code: 0.111111 (= (N-K)/K)'
b'--------------------------------'
b'infile:  ./data/data_drawnseg.txt'
b'outfile: ./data/data_rec.zip'
b'numblocks: 1'
b'assume text file as input format'
b'start decode..'
b'number of reads: 35754'
b'inner code: 0.00523018 errors on average corrected per sequence'
b'132 many erasures in block 0 of length 16383'
b'remove random sequence'
b'decoding block 0'
b'\tpart 0/6'
b'\tpart 1/6'
b'\tpart 2/6'
b'\tpart 3/6'
b'\tpart 4/6'
b'\tpart 5/6'
b'\touter code: 132 errasures, 0 errors corrected'
b'Wrote file of length 99103 Bytes'
Decoding time: 524.247209072113<jupyter_text>EvaluationLet's first extract the decoded file (data_rec.zip)<jupyter_code>with zipfile.ZipFile('./data/data_rec.zip', 'r') as zipObj:
   # Extract all the contents of zip file in decoded directory
   zipObj.extractall(path = "./data/decoded")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11871678299289236, 'fraction_numerical': 0.1530383556380867, 'mean_word_length': 2.9994878361075545, 'pattern_counts': {'":': 1, '<': 34, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '486009', 'license': 'MIT', 'n_tokens_mistral': 6947, 'n_tokens_neox': 5447, 'n_words': 1401, 'path': 'LSH_clustering.ipynb', 'repo_name': 'mersad95zd/DNA_Project', 'size': '15617'}