IGVでreadのmappingを表示する

目的

NGS(ONT MinION, long-read sequencer)の生データから興味ある遺伝子領域のみをアセンブリし、リファレンスの種や他の近縁種と配列比較することでSNPなどの検出を行いたい。さらに、short-readによるpolishが行えない(=short readのデータがない)ため、最大限readをかき集めてcoverageを稼ぎたい。

方法

QCフィルター済みのreadをdatabase、興味のある遺伝子領域の配列をqueryとしてlocal BLASTを行い、得られたreadをminimap2を用いてqueryの配列にmappingすることで、興味のある遺伝子領域のみをアセンブリできると考えた。なお、NGSデータの対象となる生物とqueryに用いた配列を持つ生物における配列は大きく変化していないことを前提とする。

問題

queryとの比較で解析対象の配列にSNPが観測された。しかし、配列を見ただけでそのSNPが真のSNPなのか、long-read sequencerのエラーによるものなのか判断することはできない。例えば、queryではTTA、解析した配列ではTGAだった場合、解析した配列で見られたSNPの箇所の全readを見てA:T:G:C = 1:1:97:1であれば真のSNPであり、A:T:G:C = 1:48:50:1であれば、エラー(あるいはpaternal・maternalの違い)である可能性が高い。
そこで、これらを判断するために重ね合わせに用いたreadがどのようにmappingされているのか視覚的に知りたい。しかし、ネットには1ページにその全てがまとまっているサイトがなかったり、ソフトのバージョンが古いためエラーが発生したりしたため、でやり方にかなり戸惑った。このページでは2020年月現在でのやり方をまとめる。

環境

macOS Catalina 10.15.5

0.現在手元にあるデータ

アセンブリまで終えると以下のファイルが揃っているはず。
・local BLASTでqueryに使用した興味のある遺伝子領域のfastaファイル
・local BLASTで回収した興味のある遺伝子領域の配列を含むreadのfastaファイル
・minimap2で作成したreadのsamファイル

1.IGVのダウンロード

Integrative Genomics Viewer (IGV)と呼ばれる有名なソフト。ChIP-seqの結果の可視化もできる。まずこのソフトをダウンロードする。
https://software.broadinstitute.org/software/igv/download

2.samtoolsのインストール

samファイルをIGVに読み込んでもらうために変換を行う。変換にはsamtoolを使用。
自分はminiconda環境があるのでそこからインストールを行う。

samtoolsのインストール

conda install -c bioconda samtools

上記の方法もしくは以下のリンク先を参考にインストールを行う。
http://kazumaxneo.hatenablog.com/entry/2019/01/04/220013

3.ファイルの準備

samtoolsを使ってsamファイルの変換を行う。
以下ターミナルから操作。

ファイルの変換

# bamファイルの作成
samtools view -bS gene.sam > gene.bam

# sorted.bamファイルの作成
samtools sort gene.bam > gene.sorted.bam

# sorted.bam.baiファイルの作成
samtools index gene.sorted.bam > gene.sorted.bam.bai

4.IGVで開く

IGVで実際のreadのmappingの様子を見る。
IGVを開き、Genomes > Load Genome from Fileでqueryに用いた配列、Files > Load from Fileでsorted.bamファイルを開く。queryの配列は一度開くとソフトに登録される。

5.IGVの見方、操作

〜編集中〜

参考サイト