macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

RNA seqのバリアントコールにも対応したABRA2

 

 次世代シーケンス(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でテストした。

Github

#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.       

                                          *|chrEBV|CMV|HBV|HCV.*|HIV.          

                                          *|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.

 

RNA

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

 

関連