シーケンシング技術の開発により、さまざまな種のゲノムを取得することが容易になっている。 NCBIゲノムデータベース(https://www.ncbi.nlm.nih.gov/genome/browse#!/)では、最大で4963の真核生物、125,679の原核生物、12,952のウイルス、10,916のプラスミド、および10,965のオルガネラゲノム配列が利用可能である; 2017年12月5日にアクセス[ref.1]。シーケンシングのエラー率は、ヒトゲノムで約0.01%である[ref.2]。しかし、ゲノム配列の品質は、特に次世代のシーケンシングプラットフォームを使用して、その後の努力によって改善されたとしても、使用されるシークエンシングプラットフォームの違いなど、さまざまな要因によって大きく異なる。さらに、Roche / 454のホモポリマーやSolexaの塩基置換など、一部のアセンブリには、使用するシーケンスプラットフォームに起因する明らかなシーケンスエラーがある[ref.3]。さらに、より多くのゲノムプロジェクトが、異なる栽培品種またはclosely relatedな種の1つのリファレンスアセンブリといくつかのリシークエンシングデータをリリースした[ref.4、5]。(一部略)
したがって、異なるアセンブリを維持および利用するための、既知のアセンブリに基づいたゲノムのアップグレード、アセンブリ、およびアノテーション付けの一般的かつ大きな要求がある。残念ながら、既知のリファレンスゲノムに基づいてゲノムのアセンブリとアノテーション転送の両方を実現する使いやすい統合ツールはほとんどない。 SAMtools / BCFtoolsやGATKなど、コンセンサスシーケンスを作成するモジュールを含むいくつかのツールはあるが、それらはいずれも各バリアントの真の対立遺伝子頻度を考慮しない。これは、偽陽性率を減らすために重要である[ref.6、7、8、9]。別のツール、Rapid Annotation Transfer Tool(RATT)はアノテーション転送に使用できるが、repetitive領域の精度は比較的低い[ref.10]。iCORNはシーケンスエラーの修正には使用できるが、アノテーションのアップグレードには使用できない[ref.11]。 Webベースのプラットフォーム(UCSC(https://genome.ucsc.edu/cgi-bin/hgLiftOver)およびGalaxy(http://usegalaxy.org))では、liftOverユーティリティを使用して異なるゲノムアセンブリバージョン間で座標を変換できるが、データベースに存在する106のゲノムのみ利用できる。遺伝子レベルでの亜種と品種間のゲノム比較の需要が高まってる。したがって、比較ゲノム解析のために、リファレンスベースのゲノムアセンブリとアノテーション転送の両方を実現することが不可欠である。残念ながら、両方の機能を実行するための不可欠なツールはほとんどない。
本研究では、ゲノムアセンブリとアノテーションプロセスで発生する問題を解決するための、リファレンスベースのゲノムアセンブリとアノテーションツールであるRGAATついて報告する。これらの問題は非常に一般的だが、Biostars(https://www.biostars.org/)とSEQanswers(http://seqanswers.com/)の2つの一般的なフォーラムを検索しても包括的なソリューションは見つからなかった。 RGAATはPerlで実装され、https://sourceforge.net/projects/rgaat/およびhttps://github.com/wushyer/RGAAT_v2からユーザーが自由に利用できる。RGAATはゲノム配列(FASTA形式)、アノテーション(GTF、GFF、GFF3、およびBED形式)、シークエンシスアライメント(SAM / BAM形式)、バリアント(VCF形式またはタブ)などのマッピングベースの新しいアセンブリ機能の入力を受け入れる。
Workflow of variant identification based on sequence alignment (SAM/BAM) in RGAAT. 論文より転載。
インストール
ubuntu18.04でテストした(docker使用、ホストmacos10.14)。
依存
cpanm Text::NSP::Measures::2D::Fisher::left
cpanm Statistics::Multtest
cpanm SVG.pm
#cpanmがなければ入れる (linux)
cd /usr/local/bin/
curl -LOk http://xrl.us/cpanm
chmod +x cpanm
git clone https://github.com/wushyer/RGAAT.git
cd RGAAT/
> perl RGAAT.pl
# perl RGAAT.pl
Author: Wanfei Liu & Shuangyang Wu
Email: <liuwf@big.ac.cn> <wushy@big.ac.cn>
Date: May 16, 2016
Version: 1.0 version
Introduction
This program can assemble and/or annotate genome for new genome and known genome upgrade using sequence alignment file (SAM or BAM format), sequence variant file (VCF format or five coloum table (tab-delimited, including chromosome, position, id, reference allele and alternative allele)) or new genome sequence file (FASTA format) based on reference genome sequence file (FASTA format) and/or annotation file (TBL, GTF, GFF, GFF3 or BED format).
Usage: perl RGAAT.pl -g genome_sequence(FASTA) -a annotation(TBL, GTF, GFF, GFF3 or BED) [-s SAM_file | -b BAM_file | -v VCF_file | -n new_genome_file] -o prefix_of_output_file.
-g: the absolute path of genome sequence file (FASTA format, required).
-a: the absolute path of genome annotation file (TBL, GTF, GFF, GFF3 or BED format, required for annotation).
-s: the absolute path of SAM file (sorted according to coordinate, it is mututally exclusive with -b, -v and -n).
-b: the absolute path of BAM file (sorted according to coordinate, it is mututally exclusive with -s, -v and -n).
-v: the absolute path of sequence variant file (VCF format or five coloum table (tab-delimited, including chromosome, position, id, reference allele and alternative allele), it is mututally exclusive with -s, -b and -n).
-n: the absolute path of new genome sequence file (FASTA format, it is mututally exclusive with -s, -b and -v).
-m: the minIdentity for BLAT comparative between reference and new genome (default value is 90).
-q: quality threshold for reads (default value is 30).
-l: read length threshold (default value is 30).
-d: read depth threshold for sequence variant (default value is 3).
-c: allele depth threshold for sequence variant (default value is 3).
-p: allele proportion threshold for sequence variant (default value is 0.5).
-fu: filter unpaired reads (yes or no, default value is no).
-fm: filter multiple mapping reads (yes or no, default value is yes).
-o: the prefix of output file (required).
Note: Genome sequence, annotation and sequence variant files should be uncompressed files. The name of annotation file should be ended with suffix *.tbl, *.gtf, *.gff, *.gff2, *.gff3 or *.bed. If you want run program multiple times with different parameters, please use "-o" option to assign different prefix name for output files. To filter false positive sequence variants, we use 3 as the default minimum reads depth and 3 as the default minimum alt allele reads coverage, because most of false positive variants come from the low depth and low alt coverage records according to our previous study. Although we filtered genome comparative result (obtained by BLAT) for reference and new genome sequence, user can manually check the *.psl file for higher accuracy if possible and only remain the most possible direction results (plus or minus). Moreover, if comparing with different assembly version or different strain of same species, we can use the default minIdentity parameter (90) for BLAT, however, if comparing with different species, we recommend using 50 as minIdentity by option -m (-m 50) for BLAT. To convert BAM file to SAM file, RGAAT need samtools.
実行方法
元のリファレンス、リードをリファレンスにmappingして得たbam、そのバリアントコールのVCFファイル、出力を指定する。アノテーションもあればそれも指定する。
perl RGAAT.pl -g referene.fasta -b align.bam -v variants.vcf -o output
出力はchange部位を記したファイル、アップデートされたゲノムなどになる。
out.svg
引用
RGAAT: A Reference-based Genome Assembly and Annotation Tool for New Genomes and Upgrade of Known Genomes
Liu W, Wu S, Lin Q, Gao S, Ding F, Zhang X, Aljohi HA, Yu J, Hu S
Genomics Proteomics Bioinformatics. 2018 Oct;16(5):373-381