次世代シーケンス（NGS）は、さまざまなアプリケーションで広く使用されるツールになっている。バリアントコールは大きな関心が寄せられている領域であり、RNAへの関心も高まっている。NGSバリアントコールパイプラインの最初のステップの1つは、シーケンスされたリードをリファレンスゲノムにアラインすることである。 bwa-memやbowtie2などの広く使用されているDNAアライナー（Langmead and Salzberg、2012; Li、2013）は、多数のリードを迅速にアラインメントし、ギャップアラインメントをサポートし、indelの識別を可能にする。速度を上げるために、これらの方法はそれぞれのリードを徹底的には調整せず、場合によっては、特に indelサイズ増加によりindelを明らかにできない場合がある。従来、バリアントコーラーは、正確にマッピングされたリードに基づいて、SNVとindelを識別する。リードが正確にマッピングされていないインスタンスでは、バリアントの検出を混乱させる可能性がある。
近年、DNAのidnel検出を改善する多くの方法が開発された。場合によっては、バリアントを含むリードがおおよそ正しい場所にマッピングされていれば、バリアントを正常に識別することができる。元のABRA実装（Mose et al、2014）では、localized assemblyプロセスを使用してリードアラインメントを調整し、その結果、アラインメントのindelを明らかにし、Freebayes（Garrison and Marth、2012）、およびsomaticコーラーのStrelka（Saunders et al、2012）などのバリアントコーラーの検出を改善した。Platypus、GATK-HaplotypeCaller（GATK-HC）、Strelka2、Scalpel（DePristo et al、2011; Kim et al、2018; Narzisi et al、2014; Rimmer et al、2014）はすべてlocalized assemblyを使用して、元のリードアラインメントに含まれる場合と含まれない場合があるインデルの検出を可能にする。同様に、最近開発されたLancetおよびMutect2（Cibulskis et al、2013; Narzisi et al、2018）は、体細胞変異コールにlocalized assemblyを使用する。
RNA-Seqは、遺伝子およびアイソフォームの発現、遺伝子融合、転写スプライシング、発現変異、RNA編集の分析を可能にする診断ツールとして非常に重要であることが証明されている。 RNAのスプライスジャンクションの存在は、スプライス対応アライナーの開発を必要とした。スプライスジャンクションにまたがるリードをマッピングできるいくつかのRNA-Seqアライナーが開発された（Dobin et al、2013; Kim et al、2013、2015; Wang et al、2010; Wu et al、2016）。ただし、些細ではないバリエーションが存在すると、RNA-Seqリードが完全にアラインメントされないことがある。例えばSunらは、標準のRNA-Seqパイプラインでは、2塩基を超える長さのindelを特定するのが困難であることを実証している (Sun, 2016)。多くの場合、元々DNA用に開発されたバリアントコールパイプラインは、RNA-Seqアラインメントの構文を処理するように変更されているが、RNA-Seq用には最適化されていない。たとえば、広く使用されているGATK（DePristo et al。、2011）では、スプライスジャンクションを削除してRNA-Seqリードアラインメントを複数のアラインメントに分割し、DNAのようなリードを生成して、DNAバリアント用に開発されたツールを使用して処理する必要がある。さらに、最近開発されたTransindel（Yang et al。、2018）は、リードの初期アラインメントをbwa-mem（DNAアライナー）に依存している。
ABRA2は、RNA-Seqデータのスプライス対応の再調整と、大幅に強化された計算パフォーマンスを提供する、元のABRA実装の更新である。 ABRA2によって生成される改善されたアライメントにより、一般的に、特にindelに対して、より正確なバリアントコールが可能になる。さらに、ABRA2は元のABRAの精度を改善し、実行時間を短縮し、ヒト全ゲノムを処理できる拡張スケーラビリティを可能にする。
簡単に言えば、元のABRA実装は、領域ごとに入力BAMファイルをローカライズして処理する。各関心の領域のリード値は、deBruijnグラフを使用して抽出およびアセンブリされる。アセンブルされたコンティグは、グローバルのコンティグプールに追加される。すべての領域でアセンブリが完了すると、bwa-memを使用してすべてのコンティグがリファレンスゲノムにアラインされる。特定のコンティグのキメラアラインメントは、アラインメントが単一のバリアントとして単純に表現できるindelの存在を明確に示している場合に結合される。すべてのコンティグがアラインメントされると、bwa-memを使用してすべてのリードがすべてのコンティグにマップされる。リードがリファレンスよりもコンティグに近い場合、リファレンスコンテキストがコンティグのアラインメントに一致するようにアップデートされる。
この方法は、遺伝子パネルやwhole exomeなど、多くのターゲットDNAシーケンスで効果的であることが証明されているが、いくつかの明らかな欠点がある。コンティグの総数が大きくなると、すべてのリードをすべてのコンティグにアラインすることが非常に遅くなり、場合によってはbwa-memが最後まで実行されないことがある。これは、多くのマウス系統の場合のように、ノイズの多いサンプルおよびリファレンスゲノムから実質的に分岐するサンプルの計算の困難さと同様に、スケーラビリティの問題を引き起こす。スケーラビリティの問題のため、ヒトの全ゲノムなどの大きなサンプルは通常、最後まで実行できない。キメラbwa-memアラインメントの後処理は、単一の単純な孤立したindel存在下ではうまく機能する可能性があるが、複数のバリアントまたは近接したテクニカルアーティファクトからなるより複雑な、またはノイズの多いイベントを適切にキャプチャできない場合がある。localized assemblyは、バリアントコールで広く使用されるようになったが、計算コストが高く、適切なバリアントを識別するために必ずしも必要ではないため、計算時間が不必要に長くなる。最後に、元のABRAもbwa-memもスプライシングを考慮していないため、ツールはRNA-Seqデータに対して最適ではない。これらの問題を改善するためにABRA2を開発した。
インストール
ubuntu18.04LTSでテストした。
#bioconda (link)
conda install -c bioconda -y abra2
> abra2
# abra2
Picked up JAVA_TOOL_OPTIONS: -Xmx4G
INFO Thu Dec 12 22:51:31 JST 2019 Abra version: 2.22
INFO Thu Dec 12 22:51:31 JST 2019 Abra params: [/root/anaconda3/share/abra2-2.22-0/abra2.jar]
Missing required input SAM/BAM file(s)
Missing required output SAM/BAM filename(s)
Missing required reference
Option Description
------ -----------
--amq <Integer> Set mapq for alignments that map
equally well to reference and an
ABRA generated contig. default of
-1 disables (default: -1)
--ca Contig anchor [M_bases_at_contig_edge,
max_mismatches_near_edge] (default:
10,2)
--cl <Integer> Compression level of output bam file
(s) (default: 5)
--cons Use positional consensus sequence when
aligning high quality soft clipping
--contigs Optional file to which assembled
contigs are written
--dist <Integer> Max read move distance (default: 1000)
--gc If specified, only reprocess regions
that contain at least one contig
containing an indel or splice
(experimental)
--gkl If specified, use the GKL Intel
Deflater.
--gtf GTF file defining exons and transcripts
--in Required list of input sam or bam file
(s) separated by comma
--in-vcf VCF containing known (or suspected)
variant sites. Very large files
should be avoided.
--index Enable BAM index generation when
outputting sorted alignments (may
require additonal memory)
--junctions Splice junctions definition file
--keep-tmp Do not delete the temporary directory
--kmer Optional assembly kmer size(delimit
with commas if multiple sizes
specified)
--log Logging level (trace,debug,info,warn,
error) (default: info)
--mac <Integer> Max assembled contigs (default: 64)
--mad <Integer> Regions with average depth exceeding
this value will be downsampled
(default: 1000)
--mapq [Integer] Minimum mapping quality for a read to
be used in assembly and be eligible
for realignment (default: 20)
--maxn [Integer] Maximum pre-pruned nodes in regional
assembly (default: 150000)
--mbq [Integer] Minimum base quality for inclusion in
assembly. This value is compared
against the sum of base qualities
per kmer position (default: 20)
--mcl [Integer] Assembly minimum contig length
(default: -1)
--mcr <Integer> Max number of cached reads per sample
per thread (default: 1000000)
--mer <Double> Min edge pruning ratio. Default value
is appropriate for relatively
sensitive somatic cases. May be
increased for improved speed in
germline only cases. (default: 0.01)
--mmr <Double> Max allowed mismatch rate when mapping
reads back to contigs (default: 0.05)
--mnf <Integer> Assembly minimum node frequency
(default: 1)
--mrn <Double> Reads with noise score exceeding this
value are not remapped.
numMismatches+(numIndels*2) <
readLength*mnr (default: 0.1)
--mrr <Integer> Regions containing more reads than
this value are not processed. Use
-1 to disable. (default: 1000000)
--msr <Integer> Max reads to keep in memory per sample
during the sort phase. When this
value is exceeded, sort spills to
disk (default: 1000000)
--no-edge-ci If specified, do not update alignments
for reads that have a complex indel
at the read edge. i.e. Do not allow
alignments like: 90M10D10I
--no-ndn If specified, do not allow adjacent N-
D-N cigar elements
--nosort Do not attempt to sort final output
--out Required list of output sam or bam file
(s) separated by comma
--rcf <Double> Minimum read candidate fraction for
triggering assembly (default: 0.01)
--ref Genome reference location
--sa Skip assembly
--sc Soft clip contig args [max_contigs,
min_base_qual,frac_high_qual_bases,
min_soft_clip_len] (default:
16,13,80,15)
--sga Scoring used for contig alignments
(match, mismatch_penalty,
gap_open_penalty,
gap_extend_penalty) (default:
8,32,48,1)
--single Input is single end
--skip If no target specified, skip
realignment of chromosomes matching
specified regex. Skipped reads are
output without modification.
Specify none to disable. (default:
GL.*|hs37d5|chr.*random|chrUn.
*|KSHV|HTLV.*|MCV|SV40|HPV.*)
--sobs Do not use observed indels in original
alignments to generate contigs
--ssc Skip usage of soft clipped sequences
as putative contigs
--sua Do not use unmapped reads anchored by
mate to trigger assembly. These
reads are still eligible to
contribute to assembly
--target-kmers BED-like file containing target
regions with per region kmer sizes
in 4th column
--targets BED file containing target regions
--threads <Integer> Number of threads (default: 4)
--tmpdir Set the temp directory (overrides java.
io.tmpdir)
--ujac If specified, use junction permuations
as contigs (Experimental - may use
excessive memory and compute times)
--undup Unset duplicate flag
--ws Processing window size and overlap
(size,overlap) (default: 400,200)
実行方法
DNA
coordinate sortされたbam、出力のbamファイル名、リファレンスゲノムを指定する。作業ディレクトリとして/tmpなどを指定する（少なくともbamと同じ程度のディスクスペースが必要）。ターゲットのbedを指定しない場合、全ゲノム領域がリアラインメント対象になる。
abra2 --in normal.bam,tumor.bam --out normal.abra.bam,tumor.abra.bam --ref hg38.fa --threads 8 --targets targets.bed --tmpdir /your/tmpdir > abra.log
- --in-vcf VCF containing known (or suspected) variant sites. Very large files （e.g. dbSNP）should be avoided.
starでmappingして得たbamを指定する。
abra2 --in star.bam --out star.abra.bam --ref hg38.fa --junctions bam --threads 8 --gtf gencode.v26.annotation.gtf --dist 500000 --sua --tmpdir /your/tmpdir > abra2.log 2>&1
- --junctions Splice junctions definition file
- --gtf file defining exons and transcripts
引用
Improved indel detection in DNA and RNA via realignment with ABRA2
Lisle E Mose, Charles M Perou, Joel S Parker
Bioinformatics, Volume 35, Issue 17, 1 September 2019, Pages 2966–2973
関連
