macでインフォマティクス

macでインフォマティクス

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

HIVディープシーケンシングのマッピングとバリアントコールパイプライン hivmmer

2018 Bioinformatics Viruses workflow HIV

　現在、いくつかの次世代シーケンシングマシンが病原体およびウイルスの研究に使用されている（Chabria et al、2014; Quin ones-Mateu et al、2014）。過去20年間に開発された多くの次世代シーケンシングプラットフォームおよびアプローチのうち、イルミナのシークエンシング技術は、主に歩留まりの向上とコストの削減により市場を支配してきた（Goodwin et al、2016）。イルミナ技術を用いたHIVサンプルのディープシーケンシングは、ウィルスの疫学、臨床的遺伝子型判定、抗レトロウィルス薬耐性の研究で頻繁に使用されている。例えば、ディープシーケンシングは、薬物耐性変異のより敏感なアッセイを提供し得る（Brumme and Poon、2016）。現在の臨床基準であるサンガーシークエンシングは、20％以下の頻度で突然変異を検出することは不可能であるか、この数値は臨床的には関連性がある可能性がある（A vilaR uteetal、2016）。ディープシーケンシングを用いた研究やこれらの臨床基準の関心は、そのシーケンシングプロトコルのエラー率の測定である。エラーは、サンプル調製（RNAゲノムのcDNAへの逆転写およびウイルスゲノムの増幅を含む）、ライブラリー調製、シークエンシングおよびbase callingにおいて生じ得る。
エラー率測定は、下流の分析において不確実性を生じる。例えば、ゲノム増幅中に導入されたエラーは、実際の突然変異と区別することは困難である。突然変異は、プロセスの初期段階で導入され、増幅され、後の段階で高い頻度で発生する可能性があるからである。 PCR中のrecombinationは、臨床的に関連する「真の」ウイルスのrecombinationと区別することは困難である。低頻度の突然変異は、シークエンシングおよびbase callingエラーと区別することが困難であり、シーケンスアライメントで正確なシーケンスオーバーラップを正確に特定することに依存するリードアライメント、アセンブリおよびハプロタイプ再構築を混乱させる可能性がある。    Beerenwinkel et al（2012）は、RT-PCRステップ中に導入されたartifactsが、ディープシーケンシングのHIVデータにおいて、個々のハプロタイプを再構成してウイルスの多様性を正確に見積もるための最大のチャレンジであると推測している。

 近年の多くのHIV研究では、Illuminaシークエンシングエラーに対して、シークエンシングエラーと区別がつかないという理由で、グローバルな頻度閾値（通常は1％）を適用してこれ以下のバリアントを除外している。このアプローチは
シーケンシングプロトコルの典型的なエラーレートを推定する必要があり、これがバリアントコールのしきい値として使用される。
シークエンシングエラー率を評価する最も一般的なアプローチは、既知のシーケンスのリードを分析することである。リードを既知シーケンスに合わせ、アライメントのミスマッチの頻度を数える。 HIVのcontentsでは、これはHIVプラスミドと既知の配列とのmixtureをシーケンシングすることによって達成できる。この研究では、このアプローチを使って、Illumina MiSeqプラットフォームによるHIV polシーケンシングの新しいパイプラインhivmmerを導入し、それを他のバリアントコールパイプラインと比較する。既存のパイプラインは、ショートリードアライナーを使用して塩基空間でHIVのリファレンス（such as HXB2; accession K03455）にアライメントするか、またはデノボアセンブリするが、hivmmerは確率的アライナーHMMER（Eddy、2011）を使用して、アミノ酸空間におけるより敏感なアライメントを行う。

 

f:id:kazumaxneo:20181113014914p:plain

Pipeline steps. Githubより

 

 hivmmerに関するツイート


インストール

Hivmmer requires Python 3.6. On 64-bit Linux, it is also possible to install hivmmer using prebuilt packages from the kantorlab Anaconda channel.

依存

  • FASTX-Toolkit 0.0.14
  • HMMER 3.2.1
  • PEAR 0.9.11

本体 Github

#公式ではcondaの仮想環境での導入を推奨しているが、ここでは直接導入
conda install -y -c kantorlab hivmmer

#dockerイメージも用意されている
docker pull kantorlab/hivmmer

> hivmmer -h

$ hivmmer -h

usage: hivmmer --id ID --fq1 FASTQ1 --fq2 FASTQ2 --ref REFERENCE [--cpu N]

               [--region REGION] [--allow-stop-codons]

               [-h|--help] [-v|--version]

 



実行方法

１、解析対象の遺伝子のprotein.fastaの準備

 

２、マルチプルアライメント及びhmmプロファイルの作成

mafft --auto input_proteins.fa > multiple_alignment.fa

#HMMERを使い、hmmプロファイルを作成
hmmbuild proteins-mafft.hmm multiple_alignment.fa

HIVのpol（HXB2 coordinates 2253-3296）（参考）に関しては、HIVデータベース（link）を元にすでにプロファイルが作成されているので、1-2の作業はスキップできる。

 

参考

HIV polに関しては、専用コマンドで作成されている。

hivmmer-trim-reference HIV1_ALL_2016_2253-3870_DNA.fasta >HIV1_ALL_2016_2253-3870_DNA.trimmed.aa.fasta

hmmbuild HIV1_ALL_2016_2253-3870_DNA.trimmed.aa.hmm HIV1_ALL_2016_2253-3870_DNA.trimmed.aa.fasta

 

３、hivmmer実行

hivmmer --id ID --fq1 pair1.fq --fq2 pair2.fq --ref proteins-mafft.hmm --cpu 4

 

 テストラン

git clone https://github.com/kantorlab/hivmmer.git
cd /hivmmer/test/
hivmmer --id ID --fq1 5VM_1.fastq --fq2 5VM_2.fastq --ref pol.hmm

 出力（hmmsearchが繰り返されるため1と2がある）

f:id:kazumaxneo:20181113020121p:plain

> head -n 20 ID.hmmsearch2.aavf

$ head -n 20 ID.hmmsearch2.aavf

##fileformat=AAVFv1.0

##fileDate=20181112 16:58:39

##source=hivmmer-0.1.2

##reference=ID.hmmsearch2.ref.consensus.fa

##INFO=<ID=AC,Number=.,Type=String,Description="Alternate Codon">

##INFO=<ID=ACC,Number=.,Type=Float,Description="Alternate Codon Count, for each Alternate Codon, in the same order as listed.">

##INFO=<ID=ACF,Number=.,Type=Float,Description="Alternate Codon Frequency, for each Alternate Codon, in the same order as listed.">

##FILTER=<ID=cov10,Description="Alternate Count < 10">

##FILTER=<ID=freq0.01,Description="Alternate Frequency < 0.01">

#CHROM GENE POS REF ALT FILTER ALT_FREQ COVERAGE INFO

pol-consensus PR 1 P P PASS 1.0 25 AC=cct;ACC=25;ACF=1.0

pol-consensus PR 2 Q Q PASS 1.0 25 AC=cag;ACC=25;ACF=1.0

pol-consensus PR 3 I I PASS 1.0 25 AC=atc;ACC=25;ACF=1.0

pol-consensus PR 4 T T PASS 1.0 25 AC=act;ACC=25;ACF=1.0

pol-consensus PR 5 L L PASS 1.0 26 AC=ctt;ACC=26;ACF=1.0

pol-consensus PR 6 W W PASS 1.0 26 AC=tgg;ACC=26;ACF=1.0

pol-consensus PR 7 Q Q PASS 1.0 29 AC=caa,cag;ACC=17,12;ACF=0.5862068965517241,0.41379310344827586

pol-consensus PR 8 R R PASS 1.0 29 AC=cga;ACC=29;ACF=1.0

pol-consensus PR 9 P P PASS 1.0 29 AC=ccc;ACC=29;ACF=1.0

pol-consensus PR 10 L L PASS 1.0 29 AC=ctc;ACC=29;ACF=1.0

 

> cat ID.hmmsearch2.ref.consensus.fa

$ cat ID.hmmsearch2.ref.consensus.fa

>pol-consensus

PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMELPGRWKPKMIGGIGGFIKVRQYD

QILIEICGHKAIGTVLVGPTPVNIIGRNLLTQIGCTLNFPISPIETVPVKLKPGMDGPKV

KQWPLTEEKIKALVEICTEMEKEGKISKIGPENPYNTPVFAIKKKDSTKWRKLVDFRELN

KRTQDFWEVQLGIPHPAGLKKKKSVTVLDVGDAYFSVPLDKDFRKYTAFTIPSINNETPG

IRYQYNVLPQGWKGSPAIFQSSLTKILEPFRKQNPDIVIYQYLDDLYVGSDLEIGQHRTK

IEELRQHLLKWGFTTPDKKHQKEPPFLWMGYELHPDKWTVQPIVLPEKDSWTVNDIQKLV

GKLNWASQIYAGIKVRQLCKLLRGTKALTEVVPLTEEAELELAENREILKEPVHGVYYDP

SKDLIAEIQKQGQGQWTYQIYQEPFKNLKTGKYARLRGAHTNDVKQLTEAVQKIATESIV

IWGKTPKFKLPIQKETWEAWWTEYWQATWIPEWEFVNTPPLVKLWYQLEKEPIIGAETF

 

引用

Measurement error and variant-calling in deep Illumina sequencing of HIV