遺伝子変異によって機能は変わる。さあ、一塩基多型を調べていこうか！

ワンピースにてお馴染みジンベイは、愛する兄を失い心が乱れてしまっているルフィに対し、「失ったものばかり数えるな」と、一喝しました。この一言が、ルフィが立ち直るきっかけの一つなり、今後の展開につながるわけです。
しかし一方で、遺伝子比較を行うような、いわゆるバイオインフォマティシャンは、機能を失った遺伝子を探しまくってます。「塩基配列が変化して、翻訳されるアミノ酸が変化して、このアミノ酸は機能に大事だから、この遺伝子がコードするタンパク質の機能は欠落しているね」、これを見つけるのが、楽しいですよね。失ったものを探すのが仕事の一つなんだよ。ジンベイには申し訳ない気持ちでいっぱいですね。

目的

くだらない話はここまでにしましょうか。
今回は一塩基多型(SNP)の調べ方を、簡潔に説明します。あくまで一例ですので、悪しからず。

なお、今回用いるゲノムデータは原核生物（バクテリア）です。真核のものとまた違うかもしれませんが、ご了承くださいね。

インストールしておくもの

・Trimmomatic
・Fastp
・SeqKit
・Smalt
・Samtools

全てcondaでインストールできます。

SmaltとSamtoolsを同じ環境で使用するとエラーが起こります。
おそらく関連プログラムがバッティングするのでしょう。
condaでインストールするのであれば、別々の環境で使うことをお勧めします。

方法

ある菌株Aの全ゲノム配列ファイルStA_WGS.fnaに菌株BのIlluminaのシーケンスデータStB_Read1.fq.gzとStB_Read2.fq.gzをマッピングし、IGVで見ていきます。

1. シーケンスデータ前処理

Illuminaで取得した場合、アダプター処理、クオリティトリミング、ダウンサンプリングが必要です。それぞれの意義はこちらに詳しく書かれています。筆者はtrimmomaticでアダプター処理、fastpでクオリティトリミング、seqkit sampleでダウンサンプリングします。アダプター処理はアダプター配列を除去するだけですので割愛しますが、クオリティトリミングとダウンサンプリングは人によってやり方は少しずつ変わってくるらしいです。そのため、二つを紹介します。

クオリティトリミング

クオリティトリミングによって、リード中の信頼性の低い塩基配列が取り除かれます。
以下のようにfastpを実行します。

# 一例
fastp -i StB_Read1.fq.gz -I StB_Read2.fq.gz -o StB_Read1.fq.qf.gz -O StB_Read2.fq.qf.gz -l 41 -q 20 -w 8 -h StB.html -j StB.json

「-i, -I」でインプットファイル、「-o, -O」でアウトプットファイル指定します。
「-l」で指定した数値未満のリードを排除します。ここでは41 bp未満のリードを排除しました。短いリードは非特異的なマッピングにつながるためです。
「-q」でクオリティボーダを指定します。シーケンスリードの塩基配列にはそれぞれ信頼値が与えられており、本オプションで指定した値をQとすると、10^(-Q/10)以下の信頼値の配列をトリミングできます。リードの信頼値に関しては別記事を参照してください。
「-h」、「-j」でクオリティトリミングのログファイルを指定できる。

ダウンサンプリング

ダウンサンプリングでは、無作為にリード削除して、データ中の総リード数を減らします。リード数を減らすことで、マッピングの処理時間や正確性を高めることができるためです。
以下のようにseqkit sampleを実行します。

# ダウンサンプリング（Read1を例にする。なお、Read2も同様である。）
seqkit sample -q VALUE -s YOUR_PREFER_NUMBER StB_Read1.fq.qf.gz > StB_Read1_down.fq.qf
# 圧縮 
gunzip StB_Read1_down.fq.qf

「-q」で取り除く割合を指定します。例えば0.1と指定すれば、総リード数の0.1倍にデータを減少させます。
「-s」で乱数を指定します。ペアエンドを使っている時に重要なオプションです。ペアエンドのファイルをそれぞれseqkit sampleにかけますが、対応しあうリードが残っていないとダメですよね。乱数を一緒にすれば、無作為にリードを減少させられるかつ対応するリードを残せます。

減少させる割合に関して。
私は、「マッピングさせるゲノムの塩基数x60~100程度」を目指します。いわゆるカバレッジがx60~100程度ってことです。例えば、リファレンスゲノムが6Mbだったら、リードを360Mb~600Mbになるまでダウンサンプリングさせます。

しかしこれはあくまでも経験則であり、数値的な根拠があるわけではないです。多くすればシーケンスエラーの割合が多くなってSNPの偽陽性が出てしまいますし、少なすぎても然り。こちらも同じような意見が書かれています。読んでみてください。

2. マッピング

シーケンスリードの処理が終わったので、マッピングしていきましょう。

# リファレンスのインデックス化
smalt index -s -k 11 StA StA_WGS.fna

# マッピング
smalt map -f samsoft -o StB_to_StA.sam StA StB_Read1_down.fq.qf.gz StB_Read2_down.fq.qf.gz

# sam -> bamに変換と、リードをソーティング
samtools sort -o StB_to_StA_sorted.bam StB_to_StA.sam

# bamファイルをインデックス化
samtools index StB_to_StA_sorted.bam

マッピング自体はこれで完了です。あとは、マッピング結果をIGVでみて、SNPを調べていきます。
IGVの使い方は別記事をご参照あれ！

終わりに

SNPを調べるための、マッピングの一例を示しました。あくまでも一例です。使うサンプルの特徴によって、変わってくると思います。ただ、本記事が少しでもお役に立てれば幸いです！