GTFとGFFフォーマット

2019 10/15追記

2020 10/13 リンク削除

GTF（General Transfer Format)）はgeneのアノテーション専用のフォーマットと定義されている。それに対してGFF（General Feature Format）はtranscriptなどにも使えるよりジェネラルなフォーマットとなっている。この違いのため、例えばUCSC genomeではgeneアノテーションファイルはgtfフォーマットでのみダウンロード可能になっている。

GFFはバージョンにも注意する。GFF3はGFFのバージョン３を意味する。GFF2とGFF3には相違がある。GFF2はGTFと同じである。

フォーマットの説明；

GTF、GFF3いずれも９のカラムからなるが、1〜8行目はGTFとGFFで同じのため、GTFを例に1-8行目を説明する。例えば以下はUCSCのgenomeデータベースからダウンロードしたバクテリアのGTFファイルの最初の1行を表示している。

chr synePCC6_refSeq start_codon 1 3 2.000000 + . gene_id "slr0611"; transcript_id "slr0611"; 

1. リファレンス名 (chromosome1など)
2. アノテーションのソース (pfam,blast2go,interpro,est)
3. feature (gene,exon,start_codon,cds,mRNA,zinc_finger,conserved_region)
4. 1-based, inclusive start coordinate (integer > 0)
5. 1-based, inclusive end coordinate (integer > 0)
6. スコア
7. リファレンスに対するgeneの向き (+はforward,-はreverse)
8. frame (0,1,2 or .)

 Scoreはアノテーションの信頼性を表す数値で、gtfでは無評価の.になるのが普通。

 その他、Gffは#をつけることでコメント行を残すことが許されているなどの違いもある。

変換する簡単なスクリプトが公開されている。

SEQanswers - View Single Post - Tab Delimited File Editors? (GFF to GTF)

追記 

圧縮してindexをつける（問い合わせするため、IGVなどのビューアもindexを必要とする）。

#1 GFF
bgzip -c hg38.gff > hg38.gff.gz
tabix -p gff hg38.gff.gz

#2 GTF - tabixはGTFに対応していないので#を消すなどの細工をしてGFFに見せかける（biostarのスレッドより)
(grep ^"#" hg38.gtf; grep -v ^"#" hg38.gtf | sort -k1,1 -k4,4n) | bgzip > hg38.gtf.gz
#tabixでGFFにindexをつける（-p gff|bed|sam|vcf）
tabix -p gff hg38.gtf.gz

 tabixとbgzipがないならhtslibライブラリをコンパイルしてパスを通す（samtoolsのインストール参照）。

GTF => GFF3（gffreadを使う）

gffread transcripts.gtf -g genome.fa -E -o transcripts.gff3

IGV 

2021 8/8 こちらにまとまった情報があります。

https://www.biostars.org/p/368167/

一番下の表が分かりやすい

https://github.com/NBISweden/AGAT/blob/master/docs/gff_to_gtf.md

 関連

GFFの標準化、修正、属性の変更、追加、要約統計、イントロンサイズ、UTR情報、隣接遺伝子との重複の修正、フォーマット変換などなど。GTFとGFFの歴史についてもまとまっています。

Tandem duplicationの検出テスト

 最後はtandem duplicationのテスト結果についてまとめる。

検証

逆位の場合と同じようにシミュレーションデータを使って検証した（read-pairは除く）。結果だけ箇条書きする。  

• read-pair法のBreakdancerは100 bp以上のtandem duplicationを全て検出した。
• Split-read法のPindelは10bp -10 kbのtandem duplicationを全て検出し、配列予測も正確だった。
• Split-read法のBreseqは10 bpのサイズのみ100%検出した。
• アセンブルのPlatypusは10bpのサイズのみ100%検出した。 アセンブルのFermikitとSvABAは10bp-10kbのtandem duplicationを全て検出した。ただし両方ともサイズを正しく検出できていなかった。
• hybird法のScanindelは300bp以上の逆位を全く検出しなかった。 

まとめ

今回の検証では100bpのtandem duplicationの検出が限界で、リードより長くなってくるとアセンブリベースの手法の精度が損なわれる結果となった。挿入で良好な成績だったScanindelが長いサイズに対応できなかったのも、アセンブリが原因と考えられる。

反面、BreakdancerとPindelはリードより長いtandem duplicationについて高い検出率を示した。 　逆位の検出の結論と同じになるが、複雑な構造変化の検出率を上げるには、複数手法の結果の整合性を取っていく必要があると考えられる。 

Inversionの検出テスト

・検証  

リアルデータでは既知の逆位変異がなかったので、シミュレーションデータだけ使って3手法のパフォーマンスを検証した。箇条書きで記す。

• read-pair法のBreakdancerは、100 bp以上の逆位を100%検出した。
• Split-read法のPindelは10bp -10 kbの逆位を100%検出し、配列予測も正確だった。
• Split-read法のBreseqは100bp以下の挿入を100%検出したが、より大きなサイズだとsmall indelと間違えて検出した。
• アセンブルのPlatypusは10bpのサイズのみ100%検出した。
• アセンブルのFermikitとSvABAは10-100bpの逆位を100%検出した。より大きな逆位も検出したが、逆位の両端を２つのsmall indelと間違えて検出した。
• hybird法のScanindelは逆位を全く検出しなかった。

以上のことから、検出可能なINversionのサイズを下表のようにまとめた。

・まとめ

　現状、アセンブリベースの手法では長い逆位を正しく検出することができない。また、ScanIndelは逆位を全く検出しなかったので、検出に対応していないと思われる。複雑な構造変化の検出率を上げるには、BreakdancerやPindelなどの手法で可能性がある部位を見積もって、アセンブリベースの手法の結果と整合性を取っていく必要があると考えられる。

large insertionの検出テスト

欠損に続き挿入も検出できるかテストしたので報告する。

検証

deletionと同じようにシミュレーションとリアルデータ両方を使って、read-pair以外の手法を検証した。論文化がまだなので、結果だけ箇条書きする。

• read-pair法のBreakdancerはシミュレーションデータでは100 bp以上の挿入を100%検出したが、リアルデータではdefault条件でランできなかった。
• Split-read法のPindel、Breseは、シミュレーションデータ、リアルデータに関わらず100bp以下の挿入を検出した。リード長以上の欠損は全く検出しなかった（論文に記載されている通り）。
• アセンブルのPlatypus、Fermikit、そしてSvABAはシミュレーションでは好成績だったが、リアルデータでは大きくバラついた。下の表の推奨ターゲットはリアルデータの結果からまとめた。
• Scanindelは100bp以下の挿入の検出はパーフェクトだったが、長くなるについれて検出率は落ちた。この傾向はシミュレーション、リアルデータ両方で共通だった。また挿入が大きなサイズだと１つを複数回検出することがあった。
• LUMPYは挿入を検出できないので除外した。

以上のことから、検出可能な挿入サイズは以下の表のようにまとめた。

まとめ

　シミュレーションでは高感度でも、リアルデータではイマイチの手法が目立った。リアルデータでベストな成績を出したScanIndelとPindelは有望と思われるが、同時に検出精度は100%に達しなかった点も明記しておく。現状ではInsertionをできる限り漏れなく捉えるには複数手法を並行して行う必要があると思われる。

large deletionの検出テスト

　以前、構造変化を検出する4つの方法を紹介した。

４つの原理を利用したツールは論文もソースコードも探せば無償で入手できるので、興味があればインストールして自分のデータを解析することも可能である（このブログでも紹介している）。ただしペーパーでは精度の高い方法として紹介されていても、リアルデータではイマイチなことも多い。理由は、データがシーケンサーやライブラリ作成でバイアスを受けていたり、コンタミリードがあったり、リードのクオリティが片方だけ悪かったり、はたまたカバレッジがそもそも薄かったりと様々である。ツールについても、OS環境の違いやメモリ要求量が厳しくてランすらできない手法もある。

　そこでここではソフトをmacにインストールして、解析結果から傾向をまとめることにした。 インストール可能で生のパフォーマンスも見えて来れば、導入しやすくなると思う。

検証条件  

クローナルなhaploidゲノムの欠損変異を検出することを想定し、10 ~ 5000bpのサイズの欠損を入れたシミュレーションデータ (250 bp x 2、50カバレッジ)とリアルデータ (4Mのバクテリア、301bp x 2、およそ30カバレッジ) を準備した。不完全欠損も検証するため、0 ~ 5000bpのサイズの欠損を25%、または50%混ぜこんだデータを用意した。ツールは、以前紹介したPindel、Breakdancer-max、Platypus、Fermikit、SvABA、Breseq、LUMPY、Scanindelを使用した。

結果

　本来は結果の詳細を載せるべきであるが、論文化できていないので手法それぞれのテスト結果を下に箇条書きする。またその結果から手法ごとに適した欠損サイズを表にまとめた。

• read-pair法のBreakdancerはシミュレーションデータの100 bp以上の完全欠損、不完全欠損を95%以上検出した。リアルデータを使うとdefault条件でランできなかった (インサートサイズがpoor qualityになる) 。

• Split-readとread-pairを組み合わせたPindelは、シミュレーションデータ、リアルデータ両方とも300 bpまでの欠損を95%以上検出した。1 kbp以上の欠損の検出率は0%だった。不完全欠損 (50%欠損、25%欠損)のシミュレーションデータを使うと、10bp以下の欠損の検出感度が大きく低下した。

• split-readとread-depthを組み合わせたBreseqは、シミュレーションデータ、リアルデータ両方とも全ての欠損を検出した。しかしシミュレーションの不完全欠損 (50%欠損、25%欠損) の検出率はゼロ%だった。

• アセンブルのFermikit、そしてSVABAは短い欠損も長い欠損も全て検出した（完全欠損のシミュレーション）。リアルデータでは結果が安定しなかった。不完全欠損 (50%欠損、25%欠損)のシミュレーションだとサイズに相関を見せずバラついた。
• アセンブルのPlatypusはリアルデータ、シミュレーションデータの完全欠損、不完全いずれも数十bp以内の短い欠損を100%検出した。

• hybrid法のScanindelは、リアルデータ、シミュレーションの完全欠損で100%欠損を検出した。ただしシミュレーションでは大きな欠損の１部位を複数回検出する傾向があった。シミュレーションの不完全欠損 (50%欠損、25%欠損)だと、50%欠損は100％検出したが、25%欠損はおよそ50%の部位だけ検出した。

• hybrid法のLUMPYは100bp以上の欠損のみ100%検出したが、サイズのずれや冗長な検出が多かった。

以上のことから、精度よく検出可能な欠損サイズは以下のようにまとめた。 

・まとめ 

　結果は短い欠損か長い欠損のみ検出するツールと、サイズによらず検出可能なツールに別れた。また、シミュレーションでは高感度でも、リアルデータでは検出感度が悪くなる傾向が強く見られた。例えばシミュレーションで万能だったFermikitは、リアルデータのlong deltionを全く検出できなかった。一方、複数手法を混合しているScanindelとBreseqはサイズによらず高感度に検出した。各アルゴリズムの長所をうまく組みあ合わせているからと考えられる。ただし不完全欠損のシミュレーションでは、Breseqは検出感度がゼロになった。scanindelも欠損がレアだと感度が下がる傾向が見られた。また、リアルデータの長い不完全な欠損を検出できた手法は１つもなかった。

　欠損変異は4つの手法いずれでも検出は可能で、他の複雑な構造変化と比べると検出は容易と考えられる。それでも検出漏れがたくさん出たことは、万能の欠損検出ツールというのは無くて、複数のツールを使い分けるのが必須なことを示唆している。特にコピー遺伝子の欠損変異 (CNVs) や、不完全欠損などの欠損を高精度に検出するには今後さらなる手法の改善が必要なことを示している。解決にはlong read情報が必須なのかもしれない。 　

　今回使ったリアルデータは30カバレッジのややカバレッジが少ないデータである。ツールのパラメータ設定、生物種、リードのサイズ、カバレッジ、クオリティにより上記の結果は変わる可能性がある。

bamファイルの分離とマージ

複数回シーケンスしたデータを統合するため、bamファイルをmergeすることがある。

gatkのチームもこの話題を取り上げており、以下のURLで見ることができる。

https://software.broadinstitute.org/gatk/documentation/article.php?id=3060

mergeするにはsamtoolsのmergeコマンドを使う。

samtools merge merged.bam 1.bam 2.bam

merged.bamが出力される。リードグループに注意する。

bamから特定のクロモソーム情報を取り出すなら、

samtools view -b recal_reads.bam I chr19 > chr19.bam

このように行う。

fasta名がchr19になる。個人的にはスペース混じりの複雑な名前でも抽出してくれるのがありがたい。

計算リソースが潤沢な環境なら、分割してパラレルにランして解析速度を速めるような用途にも使えると思われる。

追記

elprep

sam/bam関係のツールまとめ - macでインフォマティクス

split   インプット (sam/bam/cram)をクロモソームごとに分割する。unmapリードは別出力する。

elprep split input.sam output/ --output-prefix "split-sam" --output-type sam --nr-of-threads 8 --single-end 

 ソートやindexファイルがなくても動作する。シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。 

merge  splitで分割した複数のインプット (sam/bam/cram)をマージする。

elprep merge input/ output.sam --nr-of-threads 8 

シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。bamとcram出力にも対応している。

サンプル間で共通する変異と固有の変異を抽出する

2019 11/19 コマンドエラー修正（A、Bの比較なのにCの表記がある）

2020 10/20 インストール追記

以前ショートリードからindelとSNVを検出するワークフローを紹介した。

複数サンプルがある場合、上記のような方法でVCFファイルを出力した後、サンプル間で共通のSNPs、サンプルごとに固有のSNPsなどを絞り込む必要が出てくるシチュエーションは多いと思われる。やり方は複数存在するが、今回は２つ紹介する。

9/14追記

A、gatkを使って行う

SelectVariants

https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

以下の２ステップで固有/共通変異を識別する。 

１、比較するvcfファイルを結合（２サンプル）

java -jar GenomeAnalysisTK.jar -T CombineVariants -R ref.fa -genotypeMergeOptions UNIQUIFY -V sample1.vcf -V sample2.vcf -o union.vcf

 -V: vcfファイルの指定。

genotypeMergeOptions: サンプルが同じ場合に区別するために必要。INFOカラムにsuffixがつく。詳細は公式HP参考。

(Determines how we should merge genotype records for samples shared across the ROD files (UNIQUIFY| PRIORITIZE|UNSORTED|REQUIRE_UNIQUE)

-Vを増やすことで３つ以上のサンプルに対応可。

brewでgatkを入れている人はjava -jar GenomeAnalysisTK.jarの部分をgatkだけに変えて上記コマンドを打つ。

２、共通する変異の抽出

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -o intersect.vcf

３、固有の変異（= 全変異 ー 共通する変異）

java -jar GenomeAnalysisTK.jar -T SelectVariants -R ref.fa -V union.vcf -select 'set == "Intersection";' -invertSelect -o invertselect.vcf

-invertSelect: 共通しない変異

INFOの右端にサンプル名が付いており、それでどちらのサンプルか判別可能になっている。この情報を指標に簡単なスクリプトで分離することもできる。

１のステップで結合するVCFを追加することで３つ以上のデータの共通、固有変異の抽出も可能であるし、VCFフォーマットで記述されてるデータなら、他のindel toolで検出した変異をマージすることも簡単にできると思われる。

追記; gatk、freebayes、svabaはmegeできることを確認したが、platypusやFermikitなど一部のツールではVCFのフォーマットに独自の部分があり、mergeできなかった。

B、bcftoolsを使って行う

bcftoolsのisecを使う。やり方は公式ページのisecの説明部分に書かれている。

BCFtools manual page

bcftoolsを使うには、vcfファイルが圧縮され、indexがつけられている必要がある。

準備１　bgzipで圧縮する。

bgzip -c sampleA.vcf > sampleA.vcf.gz
bgzip -c sampleB.vcf > sampleB.vcf.gz
bgzip -c sampleC.vcf > sampleC.vcf.gz

-cをつけて実行することで、元のファイルを保持したままsample.vcf.gzができる。

 準備２　インデックスファイル作成 

bcftools index sampleA.vcf.gz
bcftools index sampleB.vcf.gz
bcftools index sampleC.vcf.gz

 これで準備ができた。いかに４つの例を紹介する。

１、AとBの比較

bcftools isec sampleA.vcf.gz sampleB.vcf.gz -p output

 -p: 出力ディレクトリ（必須）

outputディレクトリにRWADMEと４つのファイルができる。仕分け条件はREADMEに書いてあるが、

・１つ目: sampleAに固有の変異（アレル）

・２つ目: sampleBに固有の変異（アレル）

・３つ目: sapleAから探したAとB共通の変異（アレル）

・４つ目: sapleBから探したAとB共通の変異（アレル）

３と４はどちらのVCFをベースに作るかの違いで、検出される部位の数は完全に同じである。VCFの由来かは４つのVCFファイルそれぞれのINFOフィールドの右端に記載されている。

２、AとBとCを比較してA特異的な変異

bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -C

-C: output positions present only in the first file but missing in the others

-Cをつけると、最初のsampleA.vcf特異的な変異だけが出力される。

３、AとBとC共通の変異

bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -n=3

-n: output positions present in this many (=), this many or more (+), this many or fewer (-), or the exact same (~) files

-nをつけることで比較するVCFファイルの数を調節する。

-n=3 → ３つのVCFファイル

-n+3 → ３つ以上のVCFファイル (３を含む）

-n-3  → ３つ以下のVCFファイル(３を含む）

上記の場合、-n=3でも-n+3でも同じ結果になる。

４、A、B、C固有の変異

bcftools isec sampleA.vcf.gz sampleB.vcf.gz sampleC.vcf.gz -p output -n-1

-n-1をつけることで、sampleA.vcf、sampleB.vcf、sampleC.vcf固有の変異が３つのファイルそれぞれに出力される。

それ以外にもたくさんオプションはある。

・出力する変異の種類を調節するには、-cオプションを使う。

　　-c all: 全変異

　　-c snps: 塩基置換のみ

　　-c indels: indelのみ

・出力ファイルを制限するには-wを使う。例えばsampleAのみの結果だけで良いなら

　　-w1

・変異の種類を細かく調節するなら-e（trueを排除）と-i（trueを残す）を使う。

　　-e'MAF[0]<0.05'　select rare variants at 5% cutoff（変異率5%以下の部位は排除。-iなら55以下だけ残す）。

　　-i'QUAL>10'　QUALが１０以上を残す。

など無数にある。他は公式ページを参照。

高度なフィルタリングでないなら、grepを使い絞り込むこともできる。

例えば、以下のようなMT（下側）にだけ検出される塩基置換を検出したいとする。

この部位のVCFファイルを確認する（青線部分）。

VCFは前もってマージしてある。

（↓実行したコマンド）

gatk -T CombineVariants -R ref.fasta -genotypeMergeOptions UNIQUIFY -V WT_SNPs.vcf -V MT_SNPs.vcf -o SNPsunion.vcf 

GT:AD:DP:GQ:PL フィールドを拡大。

GT:AD:DP:GQ:PL ./. 1/1:0,12:12:36:450,36,0

 WTでは変異が全くコールされていないため、./.になっている。MTでは塩基置換が100%おきているため0,12（REF=0、ALT=12）になっている。MTにだけおきている変異をGT:AD:DP:GQ:PL　./.の行を検索するだけで絞り込める。

ホモ(GT=1/1)

grep -n "GT:AD:DP:GQ:PL\s\./\.\s1/1" SNPs_invertselect.vcf > MT_homo.vcf

ヘテロ（GT=0/1）

grep -n "GT:AD:DP:GQ:PL\s\./\.\s0/1" SNPs_invertselect.vcf > MT_hetero.vcf

VCFtoolsのvcf-isecでも同等のことができます。VCFtoolsはベン図を書くために共通/非共通の数を知りたい場合などにも活用可能です。

追記

10Mb以下のスモールゲノムならbreseqを使うのも手です。

引用 

How to compare two vcf files of two cultivars? — GATK-Forum

GATK4のインストール

#bioconda gatk4(link) gatk4-sparkもある
conda install -c bioconda gatk4

#初回はgatk-registerを走らせて実行ファイルのパスを叩く
gatk-register <path>/<to>/<your>/GenomeAnalysisTK.jar