Zenn
🧬

【バイオインフォマティクス】ヒトT2TゲノムのResequence~Annotation

2024/05/07に公開
Python
バイオインフォマティクス
tech

この記事の要約

  • 生命科学と情報学のハイブリット学問のことを「バイオインフォマティクス」という。
  • 最もメジャーな解析にResequence（≒ゲノム変異解析）がある。
  • Resequenceは間違い探しのようなもので、基準となるゲノムAと比較するゲノムBの違いを明らかにする。
  • 変異が検出できたら、その変異はどんな遺伝子に影響を及ぼしているのかAnnotationをつける。
  • 最近人間の基準ゲノムが完璧なものに進化した。
  • 従来通り、OSSベースの解析ツールを用いて一般的なResequence解析をすると変異解析自体は問題なくできる。
  • AnnotationはAnnotation用のデータのフォーマットが従来と少し変わっており、少し工夫が必要だった。
  • 当社の解析技術共有プラットフォーム「ANCAT」で無料で解析できます。
    https://www.ancatbeta.anplat.co.jp
  • Annotationデータをよしなに整形し直してくれるpythonスクリプトを公開します。
    https://github.com/ANPLAT-Co-Ltd/modify_gff_for_human_T2T

はじめに

皆さんこんにちは。
株式会社アンプラットです。
割と最近、ヒトゲノムが完全解読されました。
ざっくり説明すると、いままでは解読がどうしても難しい領域があり、人間の完璧なゲノムはわかっておらず、研究者たちは、基準ゲノムに若干の曖昧さを残したまま、ゲノム間違い探し（遺伝子多型・変異解析）を行っていました。しかし、最近染色体の端から端まで完全解読されるという教科書にのるレベルの革新がおきました。
このあたりの解説については神がかったクオリティのレビュー論文（日本語）があるので、是非読んでみてください。
https://www.jstage.jst.go.jp/article/jsbibr/3/1/3_jsbibr.2022.primer2/_html/-char/ja

この記事では、当社のインターン生「Mさん」が、既存解析手法でヒト完全ゲノム（T2T）の解析が行えるのか調べてみた結果を記します。
※ 記事は当社社員、役員でもテストランを含む入念な確認をしており、ある程度の信頼性は担保できていますが、完全ではないので、そこはご認識の上、ご覧くださいませ。

ヒトT2TゲノムのResequence~Annotation

実際にT2Tのreference配列を用いてreseq~annotationまで行った際の流れを記述する。

目次

用意したファイル

次のような形でファイルを用意する。

/data
├── raw_fastq
│   ├── SRR5204297_1.fastq.gz
│   └── SRR5204297_2.fastq.gz
├── T2T-CHM13v2.0
│   └── genes.gff.gz
└── genomes
　　  └── T2T-CHM13v2.0.fa.gz

  • fastqファイル(SRR5204297_1.fastq.gz, SRR5204297_2.fastq.gz)

    $cd /home/data/raw_fastq
#SRA Toolkitのfastq-dumpが使えるようになっている必要あり
$fastq-dump SRR5204297 --split-files
$gzip SRR5204297*

  • fastaファイル(reference): ncbi_dataset/data/GCF_009914755.1_T2T-CHM13v2.0_genomic.fna
    NCBIから入手可能(https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009914755.1/)。

    #ダウンロード後
#ファイル圧縮
$bgzip ncbi_dataset/data/GCF_009914755.1_T2T-CHM13v2.0_genomic.fna
#ファイル移動
$mv ncbi_dataset/data/GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz /home/data/genomes/T2T-CHM13v2.0.fa.gz

    このディレクトリ名やファイル名はあとでsnpEffを使うときにデフォルトで指定されているものなのでここで変えておくと楽。

  • gffファイル(annotationのデータベースビルド用, 既存のデーターベースを使う場合には不要): ncbi_dataset/data/GCF_009914755.1/genomic.gff
    こちらもfastaと同様NCBIから入手可能だが、ダウンロード後少し修正が必要。

    #ダウンロード後
#ファイル圧縮
$gzip ncbi_dataset/data/GCF_009914755.1/genomic.gff.gz
#ファイル移動
$mv ncbi_dataset/data/GCF_009914755.1/genomic.gff.gz /home/data/T2T-CHM13v2.0/

    ここまではfastaと一緒。今回使ったgffファイルにおいて必要な修正は、Githubで公開してるpyファイル(modify_chrname.py, modify_gff.py)を実行することで一気に行える。具体的な修正内容についてはgffファイルの修正についての補足説明で述べる。

    #修正実行
#python3 modify_gff.py -i input_file -o output_file
$python3 /home/apps/modify_gff.py -i /home/data/T2T-CHM13v2.0/genomic.gff.gz -o /home/data/T2T-CHM13v2.0/genes.gff && gzip /home/data/T2T-CHM13v2.0/genes.gff

    *gffをbgzipで圧縮してしまうとビルドの時にデータの一部しか読まれずうまくいかないので注意
    同じくディレクトリ名やファイル名はあとでsnpEffを使うときにデフォルトで指定されているものなのでここで変えておくと楽。

    染色体番号の修正

    さらに、NCBIからダウンロードしたfastaファイルやgffファイルは染色体名がNC_060925.1のようにAccession番号で表示されていて見にくいのでchr1のような表記に変える。
    fastaファイル、gffファイルそれぞれ modify_chrname.pyを使って以下のコマンドで変えられる、

    $python modify_chrname.py -f "変更前fastaファイルパス" -t "染色体名の対応表(tsv/csv)" -of "変更後fastaファイルパス" -bgzip "bgzipのパス"

    *bgzipを使うのでbgzipのパス名を-bgzipで指定(default: `/home/tools/bgzipでコンテナ中では指定する必要なし)

    $python modify_chrname.py -g "変更前gffファイルパス" -t "染色体名の対応表(tsv/csv)" -og "変更後fastaファイルパス"

    "染色体番号の対応表(tsv/csv)"というのは変更前と変更後の染色体名を記載したtsvまたはcsvファイルのこと。1列目に変更前の染色体名、2列目に変更後の染色体名を記入する。

    NC_060925.1 Homo sapiens isolate CHM13 chromosome 1, alternate assembly T2T-CHM13v2.0	chr1
NC_060926.1 Homo sapiens isolate CHM13 chromosome 2, alternate assembly T2T-CHM13v2.0	chr2
NC_060927.1 Homo sapiens isolate CHM13 chromosome 3, alternate assembly T2T-CHM13v2.0	chr3
︙

    変更前の染色体名をfasta/gffファイルから抜き出すモードを modify_chrname.pyに入れたのでそれを使えば変更前の染色体名が1列に並んだファイルがつくれるので、2列目にそれぞれ変更後の染色体名を追加すれば対応表ができる。

    $python "modify_chrname.pyのpath" -f "fastaファイルパス" -m txt_fasta -ot "染色体名のtxtファイルパス"

    $python "modify_chrname.pyのpath" -g "gffファイルパス" -m txt_gff -ot "染色体名のtxtファイルパス"

    コマンドを組み合わせればfastaとgff両方一気に修正することもできるが必ずしもfastaとgff中の染色体名表記が一致しているわけではないので注意。その場合、対応表を別々に作る必要があるので同時にはできない。

    • fastaファイルの染色体名修正(コンテナ内での具体例)

      変更前はfastaファイルがこんな感じになっている
       >NC_060925.1 Homo sapiens isolate CHM13 chromosome 1, alternate assembly T2T-CHM13v2.0
   CACCCtaaaccctaacccctaaccctaaccctaaccctaaccctaaccctaacccctaaaccctaaccctaaccctaaccct…
      >…の行に書いてあるのが染色体(配列)の説明。これを染色体番号の名前に変えたい。自分で一から染色体名の対応表を作っても良いが、下のコマンドを使えば変更前の染色体名をファイルから抜き出すことができる。
      #染色体名抽出
$python /home/apps/modify_chrname.py -f /home/data/genomes/T2T-CHM13v2.0.fa.gz -m txt_fasta -ot /home/data/chrom_modify/fasta.txt
$vim /home/data/chrom_modify/fasta.txt

      そうすると下のようにfastaファイル中の`>`の行が抜き出され、`fasta.txt`として下のように出力されるので

      NC_060925.1 Homo sapiens isolate CHM13 chromosome 1, alternate assembly T2T-CHM13v2.0
NC_060926.1 Homo sapiens isolate CHM13 chromosome 2, alternate assembly T2T-CHM13v2.0
NC_060927.1 Homo sapiens isolate CHM13 chromosome 3, alternate assembly T2T-CHM13v2.0
︙
      これの2列目に染色体名を入れていく。
      NC_060925.1 Homo sapiens isolate CHM13 chromosome 1, alternate assembly T2T-CHM13v2.0	chr1
NC_060926.1 Homo sapiens isolate CHM13 chromosome 2, alternate assembly T2T-CHM13v2.0	chr2
NC_060927.1 Homo sapiens isolate CHM13 chromosome 3, alternate assembly T2T-CHM13v2.0	chr3
      ︙
      下の染色体名に,が入っているのでタブ区切りにしているが、csvでも可。その場合拡張子を.txtから.csvに変えること(GFFファイルの変更はcsvでやっているので具体例はそちらを参照)。 tsv/csvファイルができたら、次のコマンドでfastaファイルを修正。
      python /home/apps/modify_chrname.py -f /home/data/genomes/T2T-CHM13v2.0.fa.gz -t /home/data/chrom_modify/fasta.txt -of /home/data/genomes/modified_T2T-CHM13v2.0.fa.gz
      すると、fastaファイルの染色体情報の部分が次のように変更される。
      >chr1
CACCCtaaaccctaacccctaaccctaaccctaaccctaaccctaaccctaacccctaaaccctaaccctaaccctaaccct…
      zgrep -i '>' modified_T2T-CHM13v2.0.fa.gz(fastaを圧縮していない場合には grep -i '>' modified_T2T-CHM13v2.0.fa)を打つと、各染色体の情報が見られて染色体名が変更されたか確認できる。
      $zgrep -i  '>' modified_T2T-CHM13v2.0.fa.gz
>chr1
>chr2
>chr3
>chr4
︙

      こんな感じで変更後の染色体名が出力されればOK。
      変更が確認できたらオリジナルのファイルと入れ替えておく(どちらも取っときたい場合はやらなくてOK. その場合、以降の解析でのインプットファイルの名前を修正版のものに変更する)
      $mv modified_T2T-CHM13v2.0.fa.gz T2T-CHM13v2.0.fa.gz

    • gffファイルの染色体名修正(コンテナ内での具体例)

      ##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build T2T-CHM13v2.0
#!genome-build-accession NCBI_Assembly:GCF_009914755.1
#!annotation-date 10/02/2023
#!annotation-source NCBI RefSeq GCF_009914755.1-RS_2023_10
##sequence-region NC_060925.1 1 248387328
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606
NC_060925.1     RefSeq  region  1       248387328       .       +       .       ID=NC_060925.1:1..248387328;Dbxref=taxon:9606;Name=1;cell-line=CHM13htert;chromosome=1;gbkey=Src;genome=chromosome;isolate=CHM13;mol_type=genomic DNA;note=haploid cell line;sex=female;tissue-type=hydatidiform mole
NC_060925.1     BestRefSeq      gene    7506    138480  .       -       .       ID=gene-LOC127239154;Dbxref=GeneID:127239154;Name=LOC127239154;description=uncharacterized LOC127239154;gbkey=Gene;gene=LOC127239154;gene_biotype=lncRNA;partial=true
NC_060925.1     BestRefSeq      ncrna   7506    138480  .       -       .       ID=rna-NR_182074.1;Parent=gene-LOC127239154;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Name=NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1     BestRefSeq      exon    138321  138480  .       -       .       ID=exon-NR_182074.1-1;Parent=rna-NR_18207
︙
      基本動作はfastaと同様。インプットファイルのオプションが少し違うので注意。
      #染色体名抽出
$python /home/apps/modify_chrname.py -g /home/data/T2T-CHM13v2.0/genes.gff.gz -m txt_gff -ot ./chrom_modify/gff.txt
#染色体名抽出書き換え
$vim /home/data/chrom_modify/gff.txt

      NC_060927.1
NC_060934.1
NC_060933.1
︙
      同様に2列目に対応するクロモソーム名を追加。今回はカンマ区切り。
      NC_060927.1,chr3
NC_060934.1,chr10
NC_060933.1,chr9
︙
      カンマ区切りの場合は拡張子を.csvに変更
      $mv gff.txt gff.csv
      つくったcsvファイルを用いてGFFファイルを修正する。
      python /home/apps/modify_chrname.py -g /home/data/T2T-CHM13v2.0/genes.gff.gz -t /home/data/chrom_modify/gff.csv -og /home/data/T2T-CHM13v2.0/modified_genes.gff.gz

      ##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build T2T-CHM13v2.0
#!genome-build-accession NCBI_Assembly:GCF_009914755.1
#!annotation-date 10/02/2023
#!annotation-source NCBI RefSeq GCF_009914755.1-RS_2023_10
##sequence-region NC_060925.1 1 248387328
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606
chr1    RefSeq  region  1       248387328       .       +       .       ID=NC_060925.1:1..248387328;Dbxref=taxon:9606;Name=1;cell-line=CHM13htert;chromosome=1;gbkey=Src;genome=chromosome;isolate=CHM13;mol_type=genomic DNA;note=haploid cell line;sex=female;tissue-type=hydatidiform mole
chr1    BestRefSeq      gene    7506    138480  .       -       .       ID=gene-LOC127239154;Dbxref=GeneID:127239154;Name=LOC127239154;description=uncharacterized LOC127239154;gbkey=Gene;gene=LOC127239154;gene_biotype=lncRNA;partial=true
chr1    BestRefSeq      ncrna   7506    138480  .       -       .       ID=rna-NR_182074.1;Parent=gene-LOC127239154;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Name=NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
chr1    BestRefSeq      exon    138321  138480  .       -       .       ID=exon-NR_182074.1-1;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels com
      一番左のクロモソームのIDの列がクロモソーム番号に変わっている。
      同様にオリジナルのファイルと入れ替えておく。
      $cd /home/data/Tw
$mv  /home/data/T2T-CHM13v2.0/modified_genes.gff.gz /home/data/T2T-CHM13v2.0/genes.gff.gz

用意したツール

  • cutadapt (fastqファイルのトリミング、クオリティコントロール(QC)用
    公式ドキュメント: https://cutadapt.readthedocs.io/en/stable/installation.html

    #install
$pip install  cutadapt

  • bwa (リードのreferenceへのマッピング)
    Github: https://github.com/lh3/bwa

    #install
$cd tools
$git clone https://github.com/lh3/bwa.git
$cd bwa
$make

    最初はbwa-mem2を用いてマッピングを行おうとしたが、こちらはbwaよりも無駄にメモリを使うらしく、indexを付ける際にメモリ不足が発生したのでbwaに変更した。

  • samtools (SAM→BAMファイルへの変換)
    Github: https://github.com/samtools/samtools

    $wget https://github.com/samtools/samtools/releases/download/1.15.1/samtools-1.15.1.tar.bz2
$tar -xvf samtools-1.15.1.tar.bz2
$cd samtools-1.15.1/
$./configure
#インストールが必要なライブラリがあれば、ここでエラーメッセージとしてでてくるので apt install "ライブラリ名"でインストールする。
$make -j12
$make install

    *一番最新のsamtools-1.19.2でやろうとしたらsortのところでエラーが出てマッピングができなかったのでver.1.15.1で実行している。

    #エラーメッセージ例
checking for gawk... no
checking for mawk... mawk
checking for gcc... gcc
checking whether the C compiler works... yes
checking for C compiler default output file name... a.out
checking for suffix of executables... 
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether the compiler supports GNU C... yes
checking whether gcc accepts -g... yes
checking for gcc option to enable C11 features... none needed
checking for grep that handles long lines and -e... /usr/bin/grep
checking for C compiler warning flags... -Wall
checking for gcc option to enable large file support... none needed
checking location of HTSlib source tree... htslib-1.19.1
checking for NcursesW wide-character library... no
checking for Ncurses library... no
checking for Curses library... no
configure: error: curses development files not found

The 'samtools tview' command uses the curses text user interface library.
Building samtools with tview requires curses/ncurses/etc development files
to be installed on the build machine; you may need to ensure a package such
as libncurses5-dev (on Debian or Ubuntu Linux) or ncurses-devel (on RPM-based
Linux distributions) is installed.

FAILED.  Either configure --without-curses or resolve this error to build
samtools successfully.

    こういうのが出てくるので下から２番目のブロックに書いてあるライブラリ的なのをインストールする。

    # コマンド例
$apt install libncurses5-dev

    再び./configureを実行するとまた足りないライブラリがあれば同様に教えてくれるので、成功するまでapt install 〇〇→./configureをくりかえす。

    さらに、fastaファイルの修正にbgzipも使いたいのでsamtoolsの中に入っているhtslibもビルドし、さらにbgzipのシンボリックリンクを使いやすいところに置いておく。

    $cd /home/tools/samtools-1.15.1/htslib-1.15.1/
$./configure
$make
$ln -s /home/tools/samtools-1.15.1/htslib-1.15.1/bgzip /home/tools/bgzip

  • gatk (変異解析)
    Github: https://github.com/broadinstitute/gatk?tab=readme-ov-file#downloading
    インストールにjavaが必要なので入っていない場合はインストールする。

    #java install
$apt install openjdk-17-jdk
#gatk install
$wget https://github.com/broadinstitute/gatk/releases/download/4.5.0.0/gatk-4.5.0.0.zip
$unzip /home/tools/gatk-4.5.0.0.zip
$rm /home/tools/gatk-4.5.0.0.zip

  • snpEff (アノテーション)
    公式ドキュメント: https://pcingola.github.io/SnpEff/snpeff/introduction/

    $wget https://sourceforge.net/projects/snpeff/files/snpEff_v4_1k_core.zip
$unzip snpEff_v4_1k_core.zip
$rm snpEff_v4_1k_core.zip

    最新版だとなぜかgff修正用のpyファイルを実行しても一部うまく読み込めないのかDBビルドがあまり上手く行かなかったので代わりにv4.1を使った。

    これらのツールが全て組み込まれたDockerfileをつくったので、それを使ってコンテナを作ればその中で全ての処理ができる。

    #まず、"reseq_docker"というファイルにDockerfileを入れる。
$mkdir reseq_docker
$cd reseq_docker
$mv /path/to/Dockerfile ./
#イメージビルド
$docker image build -t reseq .
#コンテナ作成+データディレクトリのマウント
$docker run -v /home/murakosom/reseq_T2T/Docker/mount_data:/home/data/ -itd --name reseq reseq bash
#コンテナにアクセス
$docker exec -it reseq bash

解析手順

実際にDockerコンテナの中でreseqを実行してみた手順を示す。

  • fastqファイルのQC(cutadapt)

    • アダプター配列のトリミング

      #コマンド
$mkdir trimmed_fastq
$cutadapt -e 0.1 -o ./trimmed_fastq/SRR5204297_trimmed_1.fastq.gz -p ./trimmed_fastq/SRR5204297_trimmed_2.fastq.gz ./raw_fastq/SRR5204297_1.fastq.gz ./raw_fastq/SRR5204297_2.fastq.gz

      -e: 許容するエラーの割合の最大値(~1)
      -o: トリミング後のリード1の出力ファイル
      -p: トリミング後のリード2の出力ファイル
      最後にinputファイル2つをリード1、リード2の順に並べる

      本当は-aとか-gとかでアダプター配列を指定する必要があるが、今回はしていないのでなにもトリミングされずに終了した。

      #出力
=== Summary ===

  Total read pairs processed:          7,243,729
  Pairs written (passing filters):     7,243,729 (100.0%)
  
  Total basepairs processed: 1,448,745,800 bp
    Read 1:   724,372,900 bp
    Read 2:   724,372,900 bp
  Total written (filtered):  1,448,745,800 bp (100.0%)
    Read 1:   724,372,900 bp
    Read 2:   724,372,900 bp

    • QC

      #コマンド
$mkdir QCed_fastq
$cutadapt -q 30 -o ./QCed_fastq/SRR5204297_qc_1.fastq.gz -p ./QCed_fastq/SRR5204297_qc_2.fastq.gz ./trimmed_fastq/SRR5204297_trimmed_1.fastq.gz ./trimmed_fastq/SRR5204297_trimmed_2.fastq.gz

      -q: この数値以下のクオリティ(default: phred quality + 33)の3末端をトリミングする。
      -o: トリミング後のリード1の出力ファイル
      -p: トリミング後のリード2の出力ファイル
      最後にinputファイル2つをリード1、リード2の順に並べる

      #出力
     === Summary ===
 
 Total read pairs processed:          7,243,729
 Pairs written (passing filters):     7,243,729 (100.0%)
 
 Total basepairs processed: 1,448,745,800 bp
   Read 1:   724,372,900 bp
   Read 2:   724,372,900 bp
 Quality-trimmed:             160,946,979 bp (11.1%)
   Read 1:    40,464,114 bp
   Read 2:   120,482,865 bp
 Total written (filtered):  1,287,798,821 bp (88.9%)
   Read 1:   683,908,786 bp
   Read 2:   603,890,035 bp

  • リードのマッピング (bwa)

    • referenceのindex化

      #コマンド
/home/tools/bwa/bwa index /home/data/genomes/T2T-CHM13v2.0.fa.gz
      →出力:

      • GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz.amb

      • GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz.ann

      • GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz.bwt

      • GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz.pac

      • GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz.sa

      • mapping

        #コマンド
$mkdir mapping
$mkdir SRR5204297
$/home/tools/bwa/bwa mem -t 2 /home/data/genomes/T2T-CHM13v2.0.fa.gz /home/data/QCed_fastq/SRR5204297_qc_1.fastq.gz /home/data/QCed_fastq/SRR5204297_qc_2.fastq.gz | /home/tools/samtools-1.15.1/samtools sort -@ 2 -T /home/data/mapping/SRR5204297 > /home/data/mapping/SRR5204297_sorted.bam

        →出力:

        • SRR5204297_sorted.bam

      • bamファイルのindexファイル作成 (samtools)

        #コマンド
$/home/tools/samtools-1.15.1/samtools index /home/data/mapping/SRR5204297_sorted.bam

        →出力:

        • SRR5204297_sorted.bam.bai

  • 変異解析 (GATK)

    • MarkDuplicates

      単一のDNA断片に由来していると考えられる重複リードをタグづけする。

      #作業用ディレクトリ作成
$ mkdir call_variant
#実行
$/home/tools/gatk-4.5.0.0/gatk MarkDuplicates -I /home/data/mapping/SRR5204297_sorted.bam -O /home/data/call_variant/SRR5204297_markdup.bam -M /home/data/call_variant/SRR5204297_metrics.txt

      →出力:

      • SRR5204297_markdup.bam
      • SRR5204297_metrics.txt

    • AddOrReplaceReadGroups

      解析に必要なReadGroup tagをつける。

      $/home/tools/gatk-4.5.0.0/gatk AddOrReplaceReadGroups -I /home/data/call_variant/SRR5204297_markdup.bam -O /home/data/call_variant/SRR5204297_plus_RG.bam --RGLB 1 --RGPL ILLUMINA --RGPU flowcell-barcode-lane --RGSM SRR5204297

      -I: インプットファイル
      -O: アウトプットファイル
      --RGLB: Read-Group library
      --RGPL: Read-Group platform (e.g. CAPILLARY, LS454, ILLUMINA, SOLID, HELICOS, IONTORRENT, ONT, PACBIO)
      --RGPU: Read-Group platform unit (eg. Illuminaの場合はflowcell-barcode-lane, Soildの場合はslide)
      --RGSM: Read-Group sample name [参照]

      → 出力: SRR5204297_plus_RG.bam

      #index化
$samtools index /home/data/call_variant/SRR5204297_plus_RG.bam

      →出力:

      • SRR5204297_plus_RG.bam.bai

    • HaplotypeCaller

      変異の抽出とマッピング

      #index化
$/home/tools/gatk-4.5.0.0/gatk CreateSequenceDictionary -R /home/data/genomes/T2T-CHM13v2.0.fa.gz -O /home/data/genomes/T2T-CHM13v2.0.dict
$samtools faidx /home/data/genomes/T2T-CHM13v2.0.fa.gz

      → 出力:

      • T2T-CHM13v2.0.dict
      • T2T-CHM13v2.0.fa.gz.fai
      • T2T-CHM13v2.0.fa.gz.gzi

      *fastaファイルはgzipじゃなくてbgzipで圧縮されていないとエラーがでるので注意

      $/home/tools/gatk-4.5.0.0/gatk HaplotypeCaller -R /home/data/genomes/T2T-CHM13v2.0.fa.gz -I /home/data/call_variant/SRR5204297_plus_RG.bam -O /home/data/call_variant/SRR5204297.vcf.gz

      → 出力: SRR5204297.vcf.gz

  • アノテーション (SnpEff)

    • データベースの用意

      あまり有名じゃないゲノム配列をreferenceに用いると、アノテーション用のデータベースが公開されていない場合があるので、そういったときは自分でデータベースをビルドする必要がある。

      <公開データベースの確認>

      公開されているデータベースがあるかどうかは以下のコマンドで確認できる。

      $snpEff databases | grep -i '使用したいゲノム名/種名'

      例えば、

      $snpEff databases | grep -i human

      を実行すると、

      GRCh38.p2.RefSeq    Human genome GRCh38 using RefSeq transcripts    http://downloads.sourceforge.net/project/snpeff/databases/v4_1/snpEff_v4_1_GRCh38.p2.RefSeq.zip
PhumU1.22    Pediculus_humanus    ENSEMBL_BFMPP_22_87           	http://downloads.sourceforge.net/project/snpeff/databases/v4_1/snpEff_v4_1_ENSEMBL_BFMPP_22_87.zip
PhumU2.27    Pediculus_humanus    http://downloads.sourceforge.net/project/snpeff/databases/v4_1/snpEff_v4_1_PhumU2.27.zip

      が出てきてT2Tは無さそうだなとわかる。逆にある場合はsnpEff download -v GRCh38.p2.RefSeqでデータベースがダウンロードできるらしい。

      <データベースのビルド>

      まず、referenceファイルたちが以下のような構造で存在していることを確認。
      このときおそらく設定ファイルを書き換えない限りディレクトリ名やファイル名は固定されているので注意。(gffファイルは必ずgenes.gff/genes.gff.gzという名前でデータベース名(今回だとT2T-CHM13v2.0)のディレクトリの中に存在。fastaファイルはデータベース名のディレクトリの外に別にgenomesというディレクトリをつくり、その中に"データベース名".fa.gzという名前で保存。)

      次に、configファイルの書き換えを行う。

      $vim /home/tools/snpEff/snpEff.config

      書き換える場所は以下の2点
      (1) データディレクトリの変更(L17)

      data.dir = ./data/

      から

      data.dir = /home/data/

      に変更(または、自分がfastaやgff, bamファイルを入れた先の名前)。

      (2) ビルドするデータベース情報の追加(L315あたり)

      # Homo sapiens (hg38) (UCSC) using knownGenes instead of RefSeq
hg38kg.genome : Homo_sapiens (UCSC KnownGenes)
hg38kg.reference :	http://hgdownload.cse.ucsc.edu \												# Gene information from 'table/GTF' download
   				, http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/chromFa.tar.gz \		# Genome sequence
   				, http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/refMrna.fa.gz			# mRna

      こんな感じでデータベースの情報がたくさん載っているので、そこに新しく以下のような行を追加してビルドしたいゲノムの情報を書く。

      #Genome assembly T2T-CHM13v2.0
T2T-CHM13v2.0.genome : Homo_sapiens
T2T-CHM13v2.0.reference : https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009914755.1/

      ゲノムの種類によって他にも書き換えなければいけない部分があったりする(コドン表が違うなど)。(参照: http://pcingola.github.io/SnpEff/snpeff/build_db/)
      ここまで終われば、あとは以下のコマンドでビルドできる。

      $cd /home/data/
$java -jar /home/tools/snpEff/snpEff.jar build -gff3 -v T2T-CHM13v2.0

      -gff3:ビルドに用いるファイルの種類。他にも-gff2や-gtf22などがあり指定しなくても自動的に検出してくれるらしい。
      -v: データベース名の指定

ビルドが成功するとデーターベース名のディレクトリにsnpEffectPredictor.binというファイルが追加される。

  • アノテーション実行

    $cd /home/data
$mkdir annotation
$java -Xmx16g -jar /home/tools/snpEff/snpEff.jar -v T2T-CHM13v2.0 /home/data/call_variant/SRR5204297.vcf.gz > /home/data/annotation/SRR5204297_ann_.vcf.gz
    解析に当てるメモリを-Xmx16gで増やす。
    -v: データベース名の指定
    input_sample_file > output_sample_file

gffファイルの修正についての補足説明

GFF/GTFからのDBのビルドにあたってSnpEFFが想定しているGFFのデータ表記と実際のGFFファイル中の表記が違うため、ダウンロードしたGFFファイルをそのまま使うと、以下のようにうまく遺伝子名が表記されていない部分が多く見られる。
<修正前のGFFを用いたアノテーションの結果>

SnpEffによるアノテーション結果はVCFのINFOの列に|Allele|Effect|Putative_impact|Gene Name|Gene ID|Feature type|Feature ID|Transcript biotype|Rank/total|HGVS.c|HGVS.p|cDNA_position|CDS_position|Protein_position|Distance to feature|Errors, warnings|の順番で記載されるようになっている。 (参照:SnpEff-Input & output files)
つまり、|MODIFIER|などと書いてある部分の隣が遺伝子名に当たるが、ところどころGene-gene-"〇〇"のように、所々表記が不自然になっている。

また、この時DBビルドのログを確認すると、

WARNING_GENE_NOT_FOUND: Gene 'rna-NR_029639.1' (NC_060925.1:595721-595742) does not include 'rna-MIR200B-2' (NC_060925.1:595757-595778). Created new gene 'rna-NR_029639.1.2' (NC_060925.1:595757-595778). File '/home/data/T2T-CHM13v2.0/genes.gff' line 1738	'NC_060925.1	BestRefSeq	miRNA	595757	595778	.	+	.	ID=rna-MIR200B-2;Parent=rna-NR_029639.1;Dbxref=GeneID:406984,miRBase:MIMAT0000318,HGNC:HGNC:31579,MIM:612091,miRBase:MI0000342;gbkey=ncRNA;gene=MIR200B;product=hsa-miR-200b-3p'

WARNING: Cannot find gene 'Gene_gene-LOC100129381'. Created gene 'Gene_gene-LOC100129381' for transcript 'gene-LOC100129381'. File '/data/T2T-CHM13v2.0/genes.gff' line 7449	'NC_060925.1	Curated Genomic	pseudogene	1170200	1171667ID=gene-LOC100129381;Dbxref=GeneID:100129381;Name=LOC100129381;description=ribosomal protein S7 pseudogene;gbkey=Gene;gene_name=LOC100129381;gene_biotype=pseudogene;pseudo=true'

のような警告文がたくさん出てきており、ファイル中にtranscriptは存在するのにそれに該当する遺伝子の情報がないために架空の遺伝子としてGene-gene-"〇〇"が作られてしまっていることがわかる。

さらに、実際のGFFファイルとsnpEffで想定するGFFファイルのフォーマットを比較すると、ちょこちょこ表記に乖離があることがわかる。

<snpEffが想定するGFF/GTFファイルのフォーマット>

(Building databases: GTF / GFF details)

<実際のgffファイル(修正前、一部抜粋)>

︙
NC_060925.1	BestRefSeq	gene	7506	138480	.	-	.	ID=gene-LOC127239154;Dbxref=GeneID:127239154;Name=LOC127239154;description=uncharacterized LOC127239154;gbkey=Gene;gene=LOC127239154;gene_biotype=lncRNA;partial=true
NC_060925.1	BestRefSeq	lnc_RNA	7506	138480	.	-	.	ID=rna-NR_182074.1;Parent=gene-LOC127239154;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Name=NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1	BestRefSeq	exon	138321	138480	.	-	.	ID=exon-NR_182074.1-1;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1	BestRefSeq	exon	129906	130010	.	-	.	ID=exon-NR_182074.1-2;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1	BestRefSeq	exon	109651	109660	.	-	.	ID=exon-NR_182074.1-3;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1	BestRefSeq	exon	13445	13584	.	-	.	ID=exon-NR_182074.1-4;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
NC_060925.1	BestRefSeq	exon	7506	12982	.	-	.	ID=exon-NR_182074.1-5;Parent=rna-NR_182074.1;Dbxref=GeneID:127239154,GenBank:NR_182074.1;Note=The RefSeq transcript has 72 substitutions%2C 9 non-frameshifting indels compared to this genomic sequence;exception=annotated by transcript or proteomic data;gbkey=ncRNA;gene=LOC127239154;inference=similar to RNA sequence (same species):RefSeq:NR_182074.1;product=uncharacterized LOC127239154%2C transcript variant 1;transcript_id=NR_182074.1
︙
NC_060925.1	Curated Genomic	pseudogene	144134	146717	.	-	.	ID=gene-SEPTIN14P14;Dbxref=GeneID:107105354,HGNC:HGNC:51701;Name=SEPTIN14P14;description=septin 14 pseudogene 14;gbkey=Gene;gene=SEPTIN14P14;gene_biotype=pseudogene;gene_synonym=SEPT14P14;pseudo=true
NC_060925.1	Curated Genomic	exon	144134	146717	.	-	.	ID=id-SEPTIN14P14;Parent=gene-SEPTIN14P14;Dbxref=GeneID:107105354,HGNC:HGNC:51701;gbkey=exon;gene=SEPTIN14P14;pseudo=true
NC_060925.1	Curated Genomic	pseudogene	148562	152332	.	+	.	ID=gene-CICP3;Dbxref=GeneID:100132630,HGNC:HGNC:37742;Name=CICP3;description=capicua transcriptional repressor pseudogene 3;gbkey=Gene;gene=CICP3;gene_biotype=pseudogene;pseudo=true
︙
NC_060925.1	RefSeqFE	silencer	266278	266457	.	.	.	ID=id-GeneID:129660039;Dbxref=GeneID:129660039;Note=fragment chr1:772694-772873 (GRCh37/hg19 assembly coordinates);experiment=EXISTENCE:reporter gene assay evidence [ECO:0000049][PMID:32094911];function=represses an EF-1alpha promoter that drives an FKBP-Casp9 reporter gene in an ReSE screen in K562 erythroleukemia cells {active_cell/tissue: K562};gbkey=regulatory;regulatory_class=silencer
NC_060925.1	RefSeqFE	biological_region	266278	266457	.	+	.	ID=id-GeneID:129660039-2;Dbxref=GeneID:129660039;Name=biological region;gbkey=Region;standard_name=ReSE screen-validated silencer GRCh37_chr1:772694-772873
︙
NC_060925.1	BestRefSeq	gene	771235	771293	.	-	.	ID=gene-MIR6808;Dbxref=GeneID:102466740,HGNC:HGNC:50046,miRBase:MI0022653;Name=MIR6808;description=microRNA 6808;gbkey=Gene;gene=MIR6808;gene_biotype=miRNA;gene_synonym=hsa-mir-6808
NC_060925.1	BestRefSeq	primary_transcript	771235	771293	.	-	.	ID=rna-NR_106866.1;Parent=gene-MIR6808;Dbxref=GeneID:102466740,GenBank:NR_106866.1,HGNC:HGNC:50046,miRBase:MI0022653;Name=NR_106866.1;gbkey=precursor_RNA;gene=MIR6808;product=microRNA 6808;transcript_id=NR_106866.1
NC_060925.1	BestRefSeq	exon	771235	771293	.	-	.	ID=exon-NR_106866.1-1;Parent=rna-NR_106866.1;Dbxref=GeneID:102466740,GenBank:NR_106866.1,HGNC:HGNC:50046,miRBase:MI0022653;gbkey=precursor_RNA;gene=MIR6808;product=microRNA 6808;transcript_id=NR_106866.1
NC_060925.1	BestRefSeq	miRNA	771235	771255	.	-	.	ID=rna-MIR6808;Parent=rna-NR_106866.1;Dbxref=GeneID:102466740,miRBase:MIMAT0027517,HGNC:HGNC:50046,miRBase:MI0022653;gbkey=ncRNA;gene=MIR6808;product=hsa-miR-6808-3p
NC_060925.1	BestRefSeq	exon	771235	771255	.	-	.	ID=exon-MIR6808-1;Parent=rna-MIR6808;Dbxref=GeneID:102466740,miRBase:MIMAT0027517,HGNC:HGNC:50046,miRBase:MI0022653;gbkey=ncRNA;gene=MIR6808;product=hsa-miR-6808-3p
NC_060925.1	BestRefSeq	miRNA	771267	771288	.	-	.	ID=rna-MIR6808-2;Parent=rna-NR_106866.1;Dbxref=GeneID:102466740,miRBase:MIMAT0027516,HGNC:HGNC:50046,miRBase:MI0022653;gbkey=ncRNA;gene=MIR6808;product=hsa-miR-6808-5p
NC_060925.1	BestRefSeq	exon	771267	771288	.	-	.	ID=exon-MIR6808-2-1;Parent=rna-MIR6808-2;Dbxref=GeneID:102466740,miRBase:MIMAT0027516,HGNC:HGNC:50046,miRBase:MI0022653;gbkey=ncRNA;gene=MIR6808;product=hsa-miR-6808-5p
︙

具体的な問題としては以下の通り。

例1: NCBIからダウンロードしたgffにはsnpEffでは想定されていないlnc_RNAという分類が存在する。これがtranscriptとして認識されないため、該当するexonとlncRNA,そして遺伝子(正確にはlncRNAに遺伝子風の名前をつけたもの)がうまくマッチせず、Gene-exon-"exon-id"という遺伝子が勝手に作り出される。
→解決策: 分類名をlnc_RNAからsnpEffのドキュメントに記載のあるncrna(transcriptとして認識される)に変更
修正後: 遺伝子名が"lncRNA名(LOC〇〇)"となる。

例2: PseudogeneがsnpEffではtranscriptとして扱われるにも関わらず、gffファイルではgeneとして扱っているため、対応するgeneが存在しない。これによってDBビルドでGene-gene-"Pseudogene_name"という名前の遺伝子が勝手に作り出され、遺伝子名として表記される。
→解決策: gene-"Pseudogene_name"というidを持つ遺伝子を新たにgffに追加する
修正後: 遺伝子名が"Pseudogene_name"となる。

例3: RefSeqFE (NCBI RefSeq Functional Elements, DNAの非遺伝子部分の機能領域の分類.) によって分類されるbiological_region、enhancerなどもsnpEffでは対応していなさそう。ただ、これらの領域はすごく短いこともあってそもそも変異がほぼ無いためどのように扱われるかは確認できなかった。
→解決策: 全てまとめてinter_cnsに分類する。

例4: snpEffでは存在していないprimary_transcriptという分類名が存在する。これはtranscriptとして認識されないため、対応する遺伝子のtranscript名としてTranscript-gene-"transcript(RNA)"が勝手に作られFeatureID(transcriptのIDが載るカラム)として記載される。
→解決策:primary_transcriptをtranscriptに変更。
　修正後: FeatureID がrna-"rna-id"となる。

例5: miRNA(transcriptとして認識される)のParentとして遺伝子情報が記載されていない(transcriptの情報が記載されている)ことがあり、その場合対応する遺伝子が見つからずにGene-rna-"rna-id"という遺伝子が勝手に作り出される。
→解決策: attributeにParent=gene-"直前にある遺伝子名"を追加。
　修正後: 遺伝子名が"miRNA名(MIR〇〇)"となる。

これらはDBビルド前にそれぞれ解決策に書いてあるような形でgffファイルを修正しておく必要がある。
今回使用したファイルに関してはmodify_gff.pyを実行すれば問題なくアノテーションを行えるが、他のファイルの場合にはまた別の問題が発生する可能性もあるので、その場合DBビルドのログを見て必要な修正が何か確認する必要がある。

<修正したGFFを用いたアノテーションの結果>

謎表記が消えている。

さいごに

いかがでしたか。
Zennを見ている方々にバイオインフォマティクスをメインに活動されてる方はとても少ないかとは思いますが、Mさんが検証してくれた内容が誰かしらの役に立てばと思います。

※ T2TのAnnotationデータベースがSnpEff公式に見当たらなかったため、データベースビルドを行っていますが、SnpEff公式としては基本的に非推奨な手法です。Annotationデータベースが整備されている生物種の場合は、用意されたものを使いましょう。

Discussion