Covid-19を分析して、特徴的な配列を検索した。

初めに

Covid-19の遺伝子配列を分析して、特徴的な配列を検索してみました。
人を含め生物の遺伝子には特徴的な配列があり、それは繰り返し出現します。
しかし、それらの配列は変異や欠失などによって完全に一致しないことがあります。
なので、不一致を許しながら可能な限り似ている配列を検索しました。

注意

この記事は以下の注意点があります。
・何かしらの医学的・生物的な結果をもたらすものではありません。
・私の知識不足により誤っている場合があります。
・Covid-19の全配列ではなく前半部を対象としております。

データの準備

NCBIから直接コピーした文字列を、テキストファイルとして保存しました。
そこから、加工し60文字ずつ区切りました。
(この処理は私が保存したファイル形式が煩雑だった為必要になったものです。）

import pandas
df = pd.read_csv('covid_dna.txt', sep=' ', header=None)
df = df.set_index(0)
df['row_concat'] = df[1]+df[2]+df[3]+df[4]+df[5]+df[6]

text = ""
for i in range(1, len(df)):
    if i == 100:
        break
    t = df.iloc[i, 6]
    text = text + t

text = text.upper()

len(text)
>5940
text[:20]
>'GTTCTCTAAACGAACTTTAA'

n = 99
dna = [text[i:i+n] for i in range(0, len(text), n)]

len(dna)
>60
dna[:2]
>['GTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAG',
 'GACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTTCGTCCGGGTGTGACCGAA']

ライブラリとパラメータ

ライブラリとパラメータのインポートを行いました。

import random
import math
k = 8
t = 60
N = 100

Gibs Sampling

def _randomkmers(dna,k,t):
    randomkmer = []
    for i in range(t):
        z = random.choice(range(0,len(dna[i])-k+1))
        randomkmer.append(dna[i][z:z+k])
    return randomkmer

def _consensus(motifs):
    k = len(motifs[0])
    count = _count(motifs)
    consensus = ""
    for j in range(k):
        M = 0
        frequentSymbol = ""
        for symbol in "ACGT":
            if count[symbol][j] > M:
                M = count[symbol][j]
                frequentSymbol = symbol
        consensus += frequentSymbol
    return consensus

def hamming_distance(p, q):
    count = 0
    L = len(p)
    for i in range(L):
        if p[i] != q[i]:
            count += 1
    return count

def _profile(motifs):
    profile = {}
    t = len(motifs)
    k = len(motifs[0])
    countMotifs = _count(motifs)
    for symbol in "ACGT":
        profile[symbol] = []
    for x in countMotifs:
        for y in countMotifs[x]:
            z = y / float(t)*2
            profile[x].append(z)
    return profile

def _count(motifs):
    count = {}
    k = len(motifs[0])
    for symbol in "ACGT":
        count[symbol] = [] #Genero una lista para cada nucleotido en el set count
        for j in range(k):
            count[symbol].append(1) #a cada uno le pongo un 0
    t = len(motifs)
    for i in range(t):
        for j in range(k):
            symbol = motifs[i][j] #para el simbolo de esa posicion del motivo
            count[symbol][j] += 1 #sumarle un 1 al set count en ese lugar
    return count

def _motifs(profile,dna,k):

    motifs = []
    def maxProb(pattern,profile):
        num = []
        for i in range(len(pattern)):
            proflist = profile[pattern[i]][i]
            num.append(proflist)
        product = mathprod(num)
        return product
    for x in range(len(dna)):      #for i in range t
        prob = 0
        maxpattern = ''
        for i in range(len(dna[x])-k+1):
            pattern = dna[x][i:i+k]  #kmer text
            if maxProb(pattern,profile) > prob:
                prob = maxProb(pattern,profile)
                maxpattern = pattern
        motifs.append(maxpattern)
    return motifs

def _score(motifs):
    score = 0
    consensus = _consensus(motifs)
    for i in range(len(motifs)):
        score += hamming_distance(consensus,motifs[i])
    return score

def profile_most_probable_kmer(text, k, profile):
    mostProbVal = -1
    mostProbKmer = ''

    for i in range(0, 1 + len(text) - k):
        kmer = text[i:i+k]
        probKmerVal = _pr(kmer, profile)
        if probKmerVal > mostProbVal:
            mostProbVal = probKmerVal
            mostProbKmer = kmer

    return mostProbKmer

def _pr(text, profile):
    P = 1

    for i in range(len(text)):
        P = P * profile[text[i]][i]

    return P

def mostprobprof(Text,k,profile):
    prob = 0
    maxpattern = ''
    def maxProb(pattern,profile):
        num = []
        for i in range(len(pattern)):
            proflist = profile[pattern[i]][i]
            num.append(proflist)
        product = mathprod(num)
        return product
    for i in range(len(Text)-k+1):
        pattern = Text[i:i+k]
        if prob < maxProb(pattern,profile):
            prob = maxProb(pattern,profile)
            maxpattern = pattern
    return maxpattern

def mathprod(num):
    r = 1
    for i in num:
        r *= i
    return r

#FUNCIONS

def GibbsSampler(dna,k,t,N):
    motifs = _randomkmers(dna, k, t)
    bestMotifs = motifs
    for i in range(N):
        i = random.choice(range(t))
        del motifs[i]
        profile = _profile(motifs)
        motifs.insert(i,mostprobprof(dna[i],k,profile))
        if _score(motifs) < _score(bestMotifs):
            bestMotifs = motifs
    return bestMotifs

#print(GibbsSampler(dna,k,t,N))

bestScore = 1000
for start in range(20):
    bMotifs = GibbsSampler(dna,k,t,N)
    if _score(bMotifs) < bestScore:
        bestScore = _score(bMotifs)
        bestMotifs = bMotifs

print('\n'.join(bestMotifs))
>CTTTAAAA
TGTTGCAG
AGATGGAG
...(省略)...
TGTAGAAG
TTATGATG
TAAAGAAA
TTACAAAG
TTATAAAT