はじめに
DNA配列を主成分分析するツールは複数存在するが、SAM,BAM,VCFファイルなどの独自のフォーマットが入力として固定されていることがほとんどであり、DNA配列のローデータであるFASTAファイルを入力として主成分分析を行えるツールは私が調べたところ存在していない。そこで、本記事ではPythonを使用して、リファレンスゲノムにマップ済みのFASTAファイルからSNPを抽出し、主成分分析を行うプログラムを開発し、主成分分析結果のプロットまでを行う。
入力に用いるFASTAファイル
今回の主成分分析には、AmtDBからミトコンドリアDNAを一括でダウンロードし、FASTAファイルとして保存したものを使用する。AmtDBとは、古代ヨーロッパ人のものを中心に,ミトコンドリアDNAの最新データセットを提供しているデータベースで、サンプルのメタ情報を1つのCSVファイルにまとめたメタデータを提供しており、メタデータの取得が容易に可能な他、データの内容もID,日付,位置情報,部位,文化,mtDNA ハプログループなど豊富に存在する.サンプルのデータは主にFASTAファイル形式のローデータで提供されており,一括ダウンロードやフィルタリングが可能で,非常にダウンロードしやすいインターフェイスになっている
上のリンクからAmtDBのデータダウンロードページにアクセスし、右上の「Actions」の「Download FASTA」の項目内の、「All」をクリックし、FASTAファイルをダウンロードした。主成分分析の結果のプロットした際に、説明変数で色分けする際にメタデータも用いるので、下の「Download metadata」の項目「All」から、メタデータ一式をcsvファイルでダウンロードした。
また、参照配列は改訂ケンブリッジ参照配列rCRSを用いた。当配列は以下のURLで、NCBIからFASTA形式でダウンロードした。
参照配列とのマッピングにはMafftを用いて、マッピング後のファイルをFASTA形式で保存し、今回の主成分分析プログラムの入力とする。この際、使用した参照配列の情報はFASTAファイルの一番上に保存されているものとしてプログラムを作成した。
主成分分析プログラム
from Bio import SeqIO
import pandas as pd
from sklearn.decomposition import PCA
必要なモジュールをインポートする。
#参照配列とのアラインメント済みfastaファイルを読み込む。
# 参照配列の情報はfastaファイルの最上部に保存されているものとする。
fasta_file = dir + "align_mafft.fasta"
sequences = list(SeqIO.parse(fasta_file, "fasta"))
#参照配列
ref = sequences[0]
#配列の長さ
ref_seq = ref.seq
seqLen = len(ref.seq)
ignore_base = ["n","N","-","R","r","Y","y","S","s","k","K","W","w","M","m"
,"V","v","H","h","D","d","B","b"]
主成分分析に用いるFASTAファイルを、BiopythonのSeqIOモジュールを用いてリストsequencesとして読み込む。リストの各要素は、読み込んだFASTAファイルの中の個々のサンプルのDNA配列情報である。
参照配列のDNA配列情報はFASTAファイルの最上部に保存されていることを想定しているため、リストの0番目の情報を変数refとして保存し、参照配列を変数ref_seqとして保存する。
その後、参照配列の配列の長さを変数SeqLenに保存する。参照配列とのマッピングを事前に行っているため、全てのDNA配列の長さは均等になっているが、今回は代表的に参照配列のDNA配列の長さを取得している。
また、SNPを検出する際に、あらかじめ変異としてカウントしない塩基の記号をリストignore_baseとして保存した。-はマッピング時に配列の長さを合わせるために追加された文字列、その他の記号はいずれも塩基を単一に特定できないことを示すので、変異としてはカウントしないことにした。
# snpの情報を格納する辞書
snp_dict = {}
for sample in sequences[1:]:
snp_dict[sample.id] = {} # 各配列のidをキーとする
seq = str(sample.seq)
for i in range(1,seqLen): # 参照配列と比較する
if (seq[i] in ignore_base or ref_seq[i] in ignore_base):
continue
# 変異がある場合(snpは塩基の一箇所のみが変異したものを言うので、
# 変異が連続している場合はsnpとしてカウントしない)
if(i < seqLen-1):
if seq[i] != ref_seq[i] and seq[i+1] == ref_seq[i+1]:
# snpの位置と変異型（変異後の塩基名）を辞書に追加する
column_name = str(i+1) + "_" + seq[i]
snp_dict[sample.id][column_name] = 1
else:
if seq[i] != ref_seq[i] and seq[i-1] == ref_seq[i-1]:
column_name = str(i+1) + "_" + seq[i]
snp_dict[sample.id][column_name] = 1
参照配列以外の配列と参照配列のDNA配列を前から順に比較していき,SNPを検出する.
DNAの変異はsnp_dictという辞書オブジェクトに,キーを各サンプルのidとDNA配列の各インデックス,値が各サンプルのインデックスの有無という二次元的な構造で格納する.
始めのforloopでサンプルの情報を1体ずつ読み出し,次のforloopでサンプルのDNA配列の塩基を1つずつ取り出している.
SNPのみを変異として検出するために,プログラムの25,30行目でif文を用いて,変異があったDNA配列のインデックスの次の塩基同士を比較し,塩基が一致している場合のみ,変異と見なすことにした.
配列のインデックスが最後尾の場合は,次の配列が存在しないので,例外的に変異があったインデックスの前の塩基同士を比較し,塩基が一致している場合のみ,変異と見なすことにした.
SNPを検出した場合は,検出した配列のインデックスと,置換された塩基の種類(A,T,C,Gのどれか)を取得し,column_nameという変数に保存する.
その後,サンプルの名前とcolumn_nameをキー,値を変異があったことを表す数字1として,辞書に格納する.変異塩基の種類にはA,T,C,Gの4通り存在するため,それらを区別するために,置換塩基とインデックスを組み合わせてキーとしている.
# snpのデータ行列を作成する
# 辞書からデータフレームを作成する
snp_df = pd.DataFrame.from_dict(snp_dict, orient="index")
#辞書の順番を保持
order_list = snp_dict.keys()
snp_df = snp_df.reindex(order_list)
# 変異がなかった場所は0で埋める
snp_df.fillna(0, inplace=True)
全サンプルのSNPをsnp_dictに格納したら,PandasのDataFrame.from_dict関数を用いて,辞書オブジェクトであるsnp_dictをデータフレーム形式に変形してデータ行列を作成し,snp_dfという変数に格納する.
この時点では,サンプル中で変異がなかった,あるいは変異と見なさなかった箇所は,欠損値を表すNaNが代入されているので,Pandasのfillna関数を用いて,欠損値NaNを変異がなかったことを示す数字0に置換する.
上記の処理によって作成したデータフレームsnp_dfの内容は以下のようになる。
列名はインデックス+変異塩基,行名が各サンプルのidになっており,図中の1列目の1行目は19159というサンプルの73番目の塩基が,"G"に変異しているということを表し.1列目の8行目の値は,20850の73番目の塩基が,"G"には変異していないことを表す.このデータ行列を用いて,各サンプルのミトコンドリアDNAのSNPを主成分分析する.
snp_df["sample"] = snp_df.index
samplelist1 = snp_df["sample"].tolist()
snp_df2 = snp_df.drop(columns=["sample"], axis=1)
components = 5
pca = PCA(n_components=components) # PCAオブジェクトを作成する
pca.fit(snp_df2) # 主成分分析を実行する
scores = pca.transform(snp_df2) # 主成分得点を計算する
pcadata = pd.DataFrame(
scores, columns=["PC{}".format(x + 1) for x in range(components)])
プログラム中で,主成分分析を行うためにscikit-learnモジュールのPCAクラスを用いているが,このクラスのメソッドを用いて主成分分析を行うと,元データの列の名前と行の名前が保持されない.
そのため,主成分分析を行う前に,snp_dfに"sample"という列を作成し,そこにsnp_dfの行の名前,すなわちサンプルのIDのリストを格納し,それをリスト形式に変換してsamplelistという変数に保存しておく.
このリストは後ほど行うメタデータとのデータのマージに用いる.主成分分析を行うにはデータフレーム中の全てのデータの値が定量的な値である必要があるので,snp_dfからsample列のみを削除した,snp_df2というデータフレームを作成する.
snp_df2を作成した後に,PCAクラスのfit関数を用いて,主成分分析を行う.
主成分分析を行う前に,主成分分析を行う各変数が異なる単位で測定されている場合など,各変数の範囲や分散が大きく異なる場合は,データの標準化を行う必要があるが,今回のSNPデータの各変数は全て0と1の数字のみで構成されているため,標準化のプロセスは省いた.
PCAクラスで主成分分析を行う際は,始めに特定する主成分の数を指定し,PCAクラスの引数n_componentsに指定する.今回は、２次元までのプロットしか想定していないので、適当な数字である5に設定しておく.
主成分分析を実行して主成分を特定したら,transform関数を実行して各サンプルの主成分における主成分得点を計算し,結果を変数scoresに格納する.この際,データの形式は保持され,各データの値のみが主成分得点に置き換わる.
主成分得点を格納したscoresの列名を,各主成分を表すPC1,PC2,PC3...に変更し,変数pcadataに格納する.
その後,変数pcadataにsample列を追加してsamplelistの内容を格納し,各行がどのサンプルの主成分得点を表すかを特定できるようにする.
pcadata["sample"] = samplelist1
metadata = pd.read_csv(dir + "metadata_plus_china.csv")
pcadata = pd.merge(pcadata, metadata, left_on='sample', right_on="identifier",how="inner")
pcadata.to_csv(dir + "pcadata.csv")
変数pcadataを作成したら,各サンプルのメタデータを変数metadataに保存し,metadataとpcadataをマージする.pcadataにはsample列に,metadataにはidentifier列に各サンプルのIDが格納されているので,それらをキーとしてPandasのmerge関数で内部結合を行う.
内部結合を行うので,両方のサンプルに存在するIDのみがマージ後のデータに残る.そのため,pcadataに存在するサンプルのメタデータのみをマージすることができる.
マージしたデータをpcadata.csvとして保存する.
保存したpcadata.csvの内容は以下のようになる
図の6,7列目で,各サンプルのIDがsample列とidentifier列に格納されており,両列の値が等しくなっていることが確認でき,マージが正常に行われていることが確認できる.
このデータをプロットに用いることで,国,ミトコンドリアのハプログループ,サンプルの時代などの変数でデータをカテゴライズしてプロットすることができる.
主成分分析結果のプロット
import matplotlib.pyplot as plt
#日本語対応
import japanize_matplotlib
import pandas as pd
import seaborn as sns
必要なモジュールをインポートする。
pcadata = pd.read_csv(dir + "pcadata.csv")
pcadata["mt_hg"] = pcadata["mt_hg"].str.slice(0, 1)
hg_counts = pcadata['mt_hg'].value_counts()
pcadata = pcadata[pcadata['mt_hg'].isin(hg_counts[hg_counts > 60].index)]
fig,ax = plt.subplots(figsize=(10,7))
pca_x = "PC1"
pca_y = "PC2"
hue = "mt_hg"
sns.scatterplot(x=pca_x, y=pca_y,hue = hue,data=pcadata, palette="deep")
ax.set_xlabel(pca_x)
ax.set_ylabel(pca_y)
plt.title("古代人のmtDNA主成分分析")
ax.legend(fontsize=7)
plt.show()
今回はミトコンドリアハプログループでデータをカテゴライズするが,ミトコンドリアハプログループは小分類まで含めると,U3b,U4a2a,U5a1c1など,膨大な種類が存在するため,今回は大分類のみでカテゴライズを行うことにした.
ミトコンドリアDNAのハプログループについての情報はmt_hg列に格納されている.
str.slice関数を用いて,mt_hg列の値を左から1文字目の値をmt_hg列の値として更新することで,H,U,Kなどの大分類のみの値にした.
しかし,ミトコンドリアDNAのハプログループはA〜Zまで存在するため,この状態でカテゴリ別に色分けしてプロットしても,違うカテゴリでも同じ色になってしまい,カテゴリの判別が難しくなってしまう.
そのため,Pandasのvalue_counts関数を用いて,mt_hg列中の各大分類の数をカウントし,hg_countsというリストに格納した.プログラムの11行目で,Pandasのisin関数を用いて,hg_countsの中で,値が60より大きい大分類を取得し,その大分類のサンプルのみを残すフィルタリング処理を行った.
フィルタリング処理を行った後,seabornのscatterplotメソッドを用いて,第一主成分をX軸,第二主成分をY軸として,ミトコンドリアDNAのハプログループで色分けしてプロットを行った.
プロットの結果は以下のようになった.
各ハプログループごとに独立したクラスターを形成していることが確認できた。