t-SNEでDNA配列を可視化してみた[初心者][Python]

#目的
Pythonを用いた機械学習の勉強がてら、複数のDNA配列群をt-SNEを用いてグループ分けし、可視化してみる。

#方法

- Pythonのscikit-learnを用いる。
- テスト用のDNA配列として、AT-richな配列とGC-richな配列の2群を作る。
- 特徴量には、k-mer出現頻度を用いる。k-merとは長さがkのDNA配列であり、k=3であれば、AAA, AAT, AAG... となる。k=3の時、k-merの数は4の3乗で64個となる。

#テストすること

- Perplexityによりどのようにデータが変化するか？
- どのくらいAT-rich (もしくはGC-rich)であれば、はっきりと分かれるのか？
- k-merの長さ(k)はどのくらいが良いか？
- 配列の長さはどう影響するか？

#注意点

- 筆者は初心者であり、手法やデータの解釈が正しさは保証できない。あくまで参考程度に留めておいてほしい。間違いがあれば教えていただきたい。
- 特徴量のスケール変換としてどのような手法が適切であるかは今回は考慮しない。今回は`sklearn.preprocessing.StandardScaler`を用いるが、実際は、適切な前処理を検討する必要があると思う。
- 今回はt-SNE以外の次元削除は用いない。kの値が大きい時(特徴量が多い時)には、PCA等で次元削除してからt-SNEをする方が良いかもしれない。実際はこの点も要検討だろう。

#解析の流れ

```python:tsne.py
#!/usr/bin/env python3

import random
from collections import OrderedDict
from itertools import product
import numpy as np
import pandas as pd
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt


# テスト用のDNA配列のパラメータセッティング
seq_num=500 # 作るテスト配列数
seq_len=500 # テスト配列の長さ
AT_percent=53 # テスト配列のAT%
seqs=OrderedDict()

# テスト用のDNA配列を作成
orig_seq  = 'AT'* int((seq_len * AT_percent / 200))       # AT-richな配列を作成
orig_seq += 'GC'* int((seq_len * (100 - AT_percent) / 200))
orig_seq = [ i for i in orig_seq ]
for i in range(0, seq_num):
    random.shuffle(orig_seq)
    seqs['seq'+str(i)]=''.join(orig_seq)

orig_seq  = 'GC'* int((seq_len * AT_percent / 200))       # GC-richな配列を作成
orig_seq += 'AT'* int((seq_len * (100 - AT_percent) / 200))
orig_seq = [ i for i in orig_seq ]
for i in range(seq_num, seq_num*2):
    random.shuffle(orig_seq)
    seqs['seq'+str(i)]=''.join(orig_seq)

# k-merの一覧を作る
nt=['A', 'T', 'G', 'C']
kmer_list=[''.join(i) for i in product(nt, repeat=3)] # k=3

# k-mer出現頻度を計算するための関数
def bin_slide(string, bin):
    l=[]
    for i in range(0, len(string)-bin+1):
        l.append(string[i:i+bin])
    return l

def calculate_ker_freq(string, bin):
    seq=string.upper()
    seq_kmers=bin_slide(seq, bin)
    kmer_freq=[]
    for k in kmer_list:
        kmer_freq.append(seq_kmers.count(k)/len(seq_kmers))
    return kmer_freq

# テスト配列からk-mer出現頻度を計算
df=pd.DataFrame(index=kmer_list, columns=[])
for s in seqs:
    df[s]=calculate_ker_freq(seqs[s], 3) # k=3
nda=df.T.values

# 前処理として、特徴量をスケール変換
sc=StandardScaler()
nda_std=sc.fit_transform(nda)

# t-SNE
tsne=TSNE(n_components=2)
tsne_nda=tsne.fit_transform(nda_std)
    
# 作図
fig=plt.figure(figsize=(4,4))
ax=fig.add_subplot(1,1,1)
ax.scatter(x=tsne_nda[0:seq_num,0], y=tsne_nda[0:seq_num,1], c='dodgerblue', label='AT-rich', s=15, alpha=0.25, edgecolors=None)
ax.scatter(x=tsne_nda[seq_num:seq_num*2,0], y=tsne_nda[seq_num:seq_num*2,1], c='orangered', label='GC-rich', s=15, alpha=0.25, edgecolors=None)
ax.set_title('t-SNE')
ax.legend()
ax.set_xlabel('t-SNE 1')
ax.set_ylabel('t-SNE 2')
plt.subplots_adjust(left=0.15, right=0.85, bottom=0.15, top=0.85)
plt.savefig('plot_out.pdf')

```

#とりあえずやってみる
ATが53%のランダムなDNA配列(500配列)と、GCが53%のランダムなDNA配列(500配列)を可視化(上記のスクリプト)した結果。両者が明確に分離されている。すごい。
![plot_out.jpg](https://qiita-image-store.s3.amazonaws.com/0/341384/f8f9c81d-6fe0-2e6b-8ae6-59aeccdf81aa.jpeg)

#どのくらいのAT%(GC%)の違いを区別できるか？
異なるAT%(GC%)の配列を作成し、可視化してみた。
また、[こちら](https://deepage.net/machine_learning/2017/03/08/tsne.html)によると、perplexityというパラメータが重要だそうである。そこで、同時にperplexityによる違いも観察してみる。
perplexityは10から50くらいが適切(デフォルトは30)だそうだが、今回はあえて5や100といった極端な数字を使ってみた。

結果、ATが60%の配列群と40%の配列群はかなり明確に分かれる。
ATが52%の配列群と48%の配列群もある程度明確に分かれる。
当然だが、ATが50%の配列群とGCが50%の配列群は分かれない。
![plot_out_AT_per.jpg](https://qiita-image-store.s3.amazonaws.com/0/341384/54b96edc-1d48-ad5c-0ff1-cafbfea131b4.jpeg)

<details><summary>コードはこちら</summary><div>

```python
#!/usr/bin/env python3

import random
from collections import OrderedDict
from itertools import product
import numpy as np
import pandas as pd
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt


# generate test seqs; parameters
seq_num=500
seq_len=500

# defs for calculation of k-mer occurrence
def bin_slide(string, bin):
    l=[]
    for i in range(0, len(string)-bin+1):
        l.append(string[i:i+bin])
    return l

def calculate_ker_freq(string, bin):
    seq=string.upper()
    seq_kmers=bin_slide(seq, bin)
    kmer_freq=[]
    for k in kmer_list:
        kmer_freq.append(seq_kmers.count(k)/len(seq_kmers))
    return kmer_freq

# prepare k-mer list
nt=['A', 'T', 'G', 'C']
kmer_list=[''.join(i) for i in product(nt, repeat=3)] # k=3

# plt
fig=plt.figure(figsize=(12,9))

n=0
for AT_percent in [60, 52, 50]:
    
    # generate test seqs
    seqs=OrderedDict()
    orig_seq  = 'AT'* int((seq_len * AT_percent / 200))       # generate AT-rich seqs
    orig_seq += 'GC'* int((seq_len * (100 - AT_percent) / 200))
    orig_seq = [ i for i in orig_seq ]
    for i in range(0, seq_num):
        random.shuffle(orig_seq)
        seqs['seq'+str(i)]=''.join(orig_seq)

    orig_seq  = 'GC'* int((seq_len * AT_percent / 200))       # generate GC-rich seqs
    orig_seq += 'AT'* int((seq_len * (100 - AT_percent) / 200))
    orig_seq = [ i for i in orig_seq ]
    for i in range(seq_num, seq_num*2):
        random.shuffle(orig_seq)
        seqs['seq'+str(i)]=''.join(orig_seq)

    # calculation of k-mer occurrence of test seqs
    df=pd.DataFrame(index=kmer_list, columns=[])
    for s in seqs:
        df[s]=calculate_ker_freq(seqs[s], 3)
    nda=df.T.values

    # standarize features
    sc=StandardScaler()
    nda_std=sc.fit_transform(nda)

    # t-SNE parameters
    for perplexity_value in [5, 30, 50, 100]:
        
        # t-SNE
        tsne=TSNE(n_components=2, perplexity=perplexity_value)
        tsne_nda=tsne.fit_transform(nda_std)
        
        # plt
        n += 1
        ax=fig.add_subplot(3,4,n)
        ax.scatter(x=tsne_nda[0:seq_num,0], y=tsne_nda[0:seq_num,1], c='dodgerblue', label='AT-rich', s=15, alpha=0.25, edgecolors=None)
        ax.scatter(x=tsne_nda[seq_num:seq_num*2,0], y=tsne_nda[seq_num:seq_num*2,1], c='orangered', label='GC-rich', s=15, alpha=0.25, edgecolors=None)
        ax.set_title('t-SNE, AT%%=%d, perplexity=%d' % (AT_percent, perplexity_value))
        ax.legend(fontsize=5)

plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95)
plt.savefig('plot_out.pdf')
```
</div></details>

# k-merのkの値による違いは？
kの値を1から5を用いて行ってみた。
配列は全て500-nt、AT(GC)%は52%である。
2群の違いはAT(GC)%であるので、k=1であると当然2クラスに分かれる。k=2,3ではある程度明確に分かれるが、k=4,5では微妙である。
![plot_out_kmer.jpg](https://qiita-image-store.s3.amazonaws.com/0/341384/9e183482-d468-9baf-50dc-fbb248682554.jpeg)

<details><summary>コードはこちら</summary><div>

```python
#!/usr/bin/env python3

import random
from itertools import product
from collections import OrderedDict
import numpy as np
import pandas as pd
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt


# generate test seqs
seq_num=500
seq_len=500
AT_percent=52
seqs=OrderedDict()

orig_seq  = 'AT'* int((seq_len * AT_percent / 200))       # generate AT-rich seqs
orig_seq += 'GC'* int((seq_len * (100 - AT_percent) / 200))
orig_seq = [ i for i in orig_seq ]
for i in range(0, seq_num):
    random.shuffle(orig_seq)
    seqs['seq'+str(i)]=''.join(orig_seq)

orig_seq  = 'GC'* int((seq_len * AT_percent / 200))       # generate GC-rich seqs
orig_seq += 'AT'* int((seq_len * (100 - AT_percent) / 200))
orig_seq = [ i for i in orig_seq ]
for i in range(seq_num, seq_num*2):
    random.shuffle(orig_seq)
    seqs['seq'+str(i)]=''.join(orig_seq)


# prepare k-mer lists
nt=['A', 'T', 'G', 'C']
kmer_list1=nt
kmer_list2=[''.join(i) for i in product(nt, repeat=2)]
kmer_list3=[''.join(i) for i in product(nt, repeat=3)]
kmer_list4=[''.join(i) for i in product(nt, repeat=4)]
kmer_list5=[''.join(i) for i in product(nt, repeat=5)]
kmer_lists=[kmer_list1, kmer_list2, kmer_list3, kmer_list4, kmer_list5]

# defs for calculation of k-mer occurrence
def bin_slide(string, bin):
    l=[]
    for i in range(0, len(string)-bin+1):
        l.append(string[i:i+bin])
    return l

def calculate_ker_freq(string, bin):
    seq=string.upper()
    seq_kmers=bin_slide(seq, bin)
    kmer_freq=[]
    for k in kmer_list:
        kmer_freq.append(seq_kmers.count(k)/len(seq_kmers))
    return kmer_freq

# plt
fig=plt.figure(figsize=(12,15))

n=0
for k in range(5):
    
    # calculation of k-mer occurrence of test seqs
    kmer_list=kmer_lists[k]
    df=pd.DataFrame(index=kmer_list, columns=[])
    for s in seqs:
        df[s]=calculate_ker_freq(seqs[s], k+1)
    nda=df.T.values

    # standarize features
    sc=StandardScaler()
    nda_std=sc.fit_transform(nda)

    # t-SNE parameters
    for perplexity_value in [5, 30, 50, 100]:
        
        # t-SNE
        tsne=TSNE(n_components=2, perplexity=perplexity_value)
        tsne_nda=tsne.fit_transform(nda_std)
        
        # plt
        n += 1
        ax=fig.add_subplot(5,4,n)
        ax.scatter(x=tsne_nda[0:seq_num,0], y=tsne_nda[0:seq_num,1], c='dodgerblue', label='AT-rich', s=15, alpha=0.25, edgecolors=None)
        ax.scatter(x=tsne_nda[seq_num:seq_num*2,0], y=tsne_nda[seq_num:seq_num*2,1], c='orangered', label='GC-rich', s=15, alpha=0.25, edgecolors=None)
        ax.set_title('t-SNE, k=%d, perplexity=%d' % (k+1, perplexity_value))
        ax.legend(fontsize=5)

plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95)
plt.savefig('plot_out.pdf')
```
</div></details>


# 配列の長さの影響は？
k-merのkの値を決める時には、配列の長さが大きく影響することが知られている。そこで、配列の長さとkの関係を見てみた。
なお、perplexityは50、AT(GC)は52%に固定した。

やはり、配列が長くなればなるほど、大きいkでも分離可能となっている。逆に、短い配列(200-ntの配列)では、k=2でも分かれない。
![plot_out_seq_len.jpg](https://qiita-image-store.s3.amazonaws.com/0/341384/acc482f1-85eb-8997-8337-f399785b6fa3.jpeg)

<details><summary>コードはこちら</summary><div>

```python
#!/usr/bin/env python3

import random
from collections import OrderedDict
from itertools import product
import numpy as np
import pandas as pd
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt


# generate test seqs; parameters
seq_num=500
AT_percent=52

# defs for calculation of k-mer occurrence
def bin_slide(string, bin):
    l=[]
    for i in range(0, len(string)-bin+1):
        l.append(string[i:i+bin])
    return l

def calculate_ker_freq(string, bin):
    seq=string.upper()
    seq_kmers=bin_slide(seq, bin)
    kmer_freq=[]
    for k in kmer_list:
        kmer_freq.append(seq_kmers.count(k)/len(seq_kmers))
    return kmer_freq

# prepare k-mer list
nt=['A', 'T', 'G', 'C']
kmer_list2=[''.join(i) for i in product(nt, repeat=2)]
kmer_list3=[''.join(i) for i in product(nt, repeat=3)]
kmer_list4=[''.join(i) for i in product(nt, repeat=4)]
kmer_lists=[kmer_list2, kmer_list3, kmer_list4]

# plt
fig=plt.figure(figsize=(9,9))
perplexity_value=50

n=0
for seq_len in [200, 500, 2500]:
    n += 1
    u=0
    
    # generate test seqs
    seqs=OrderedDict()
    orig_seq  = 'AT'* int((seq_len * AT_percent / 200))       # generate AT-rich seqs
    orig_seq += 'GC'* int((seq_len * (100 - AT_percent) / 200))
    orig_seq = [ i for i in orig_seq ]
    for i in range(0, seq_num):
        random.shuffle(orig_seq)
        seqs['seq'+str(i)]=''.join(orig_seq)
    
    orig_seq  = 'GC'* int((seq_len * AT_percent / 200))       # generate GC-rich seqs
    orig_seq += 'AT'* int((seq_len * (100 - AT_percent) / 200))
    orig_seq = [ i for i in orig_seq ]
    for i in range(seq_num, seq_num*2):
        random.shuffle(orig_seq)
        seqs['seq'+str(i)]=''.join(orig_seq)
    
    # t-SNE parameters
    for k in [2, 3, 4]:
    
        # calculation of k-mer occurrence of test seqs
        kmer_list=kmer_lists[k-2]
        df=pd.DataFrame(index=kmer_list, columns=[])
        for s in seqs:
            df[s]=calculate_ker_freq(seqs[s], k)
        nda=df.T.values

        # standarize features
        sc=StandardScaler()
        nda_std=sc.fit_transform(nda)

        # t-SNE
        tsne=TSNE(n_components=2, perplexity=perplexity_value)
        tsne_nda=tsne.fit_transform(nda_std)
        
        # plt
        ax=fig.add_subplot(3,3,n+(3*u))
        ax.scatter(x=tsne_nda[0:seq_num,0], y=tsne_nda[0:seq_num,1], c='dodgerblue', label='AT-rich', s=15, alpha=0.25, edgecolors=None)
        ax.scatter(x=tsne_nda[seq_num:seq_num*2,0], y=tsne_nda[seq_num:seq_num*2,1], c='orangered', label='GC-rich', s=15, alpha=0.25, edgecolors=None)
        ax.set_title('t-SNE, seq_len=%d, k=%d' % (seq_len, k))
        ax.legend(fontsize=5)
        u += 1

plt.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95)
plt.savefig('plot_out.pdf')
```
</div></details>

#感想

- t-SNEは強力にDNA配列の違いを捉え、可視化することができることを再確認した。
- 今回はAT(GC)含有量が明確に異なる配列を機械的に作成し、t-SNEを行った。これらの明確な違いでさえ、条件を変えると異なる結果が得られる。
- つまり、データの質・特徴量の選択・t-SNEのパラメータ、の全てが大事であるという、当たり前のことが理解できた。
- どのような違いがあるかわからない、より複雑な実際の生物データを用いる解析では、様々なパラメータを用いて実行し、どの条件が適切であるか検討することが必要な気がした。

#環境
Ubuntu 18.04
Python 3.7.2
pandas 0.24.1
numpy 1.15.4 
matplotlib 3.0.2 
scikit-learn 0.20.2
scipy 1.1.0