macでインフォマティクス

macでインフォマティクス

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

バクテリアのIS検出ツール IS_mapper

2019 2/19 インストールの流れを修正 

2021 8/11 condaインストール追記, help更新

 

見つけたいIS配列や抗生物質耐性カセット配列をあらかじめ入力することで、ペアエンド情報を使いISの位置を検出してくれるツール。バクテリア用に設計されており、macbook airなどのlaptopでも高速に動作する。トランスポゾンやマーカー遺伝子でタギングした系の位置の検出にも力を発揮する。論文では予測結果をPCRで確認しており、false callが非常に少ないことも主張されている。

 

原理

BWA-MEMを用いて、ペアエンドのイルミナリードをISクエリ配列にマッピングする 。得られたアラインメントファイル(SAM形式)から、SAMtools view(v0.1.19以降)を用いて、ISクエリ配列の末端にペアがマッピングされていないリード(すなわち、ISの直下に位置する配列を表すリード)を抽出し、SAMフラグに基づいてリードを検索する。具体的には、フラグ「-f 36」を用いて左フランキングリード(入力配列を左から右とする)を、フラグ「-F 40 -f 4」を用いて右フランキングリードを抽出し、別々のBAMファイルに格納した後、BedToolsを用いてFASTQ形式に変換する。Samblaster (v0.1.21) [30] を用いて、SAMファイルから、ISの末端にマッピングされ、隣接する配列にまで及んでいるリード(すなわち「ソフトクリップ」リード、図1b)を抽出する。結果として得られたFASTQファイルを、BioPythonを用いてフィルターをかけ、指定したサイズ範囲(デフォルトでは5~30bp)に収まるリードのうち、ソフトクリップされた部分を抽出する。得られた配列を左右のフランキング配列に分け、それぞれをBWA-MEMを用いて参照ゲノムまたはアセンブリマッピングすることで、解析対象のゲノムにおけるクエリISの位置を特定する。

 

インストール

依存

  • Python v2.7.5
  • BioPython v1.63
  • BWA v0.7.5a
  • Samtools v0.1.19
  • Bedtools v2.20.1 - 
  • BLAST+ v2.2.28 

依存ツールはバージョンアップしているが、最新バージョンでも問題なく動作するようである。全てpipとbrewでインストール可(例: brew install samtools & pip install pysam)。

本体 Github

git clone https://github.com/jhawkey/IS_mapper/
pip install IS_mapper/

#conda
mamba install -c bioconda ismapper -y

ismap --version

$ ismap --version

ismap 2.0

ismap -h

$ ismap -h

usage: ismap [--version] --reads READS [READS ...] --queries QUERIES

             [QUERIES ...] --reference REFERENCE [REFERENCE ...]

             [--output_dir OUTPUT_DIR] [--log LOG] [--help_all HELP_ALL]

 

Basic ISMapper options:

  --version             show program's version number and exit

  --reads READS [READS ...]

                        Paired end reads for analysing (can be gzipped)

  --queries QUERIES [QUERIES ...]

                        Multifasta file for query gene(s) (eg: insertion

                        sequence) that will be mapped to.

  --reference REFERENCE [REFERENCE ...]

                        Reference genome for typing against in genbank format

  --output_dir OUTPUT_DIR

                        Location for all output files (default is current

                        directory).

  --log LOG             Prefix for log file. If not supplied, prefix will be

                        current date and time.

  --help_all HELP_ALL   Display extended help

compiled_table.py -h

$ compiled_table.py -h

usage: compiled_table.py [-h] --tables TABLES [TABLES ...] --reference

                         REFERENCE --query QUERY [--gap GAP] [--cds CDS]

                         [--trna TRNA] [--rrna RRNA] [--imprecise IMPRECISE]

                         [--unconfident UNCONFIDENT] --out_prefix OUT_PREFIX

 

Create a table of IS hits in all isolates for ISMapper

 

optional arguments:

  -h, --help            show this help message and exit

  --tables TABLES [TABLES ...]

                        tables to compile

  --reference REFERENCE

                        gbk file of reference

  --query QUERY         fasta file for insertion sequence query for

                        compilation

  --gap GAP             distance between regions to call overlapping, default

                        is 0

  --cds CDS             qualifier containing gene information (default

                        product). Also note that all CDS features MUST have a

                        locus_tag

  --trna TRNA           qualifier containing gene information (default

                        product). Also note that all tRNA features MUST have a

                        locus_tag

  --rrna RRNA           qualifier containing gene information (default

                        product). Also note that all rRNA features MUST have a

                        locus_tag

  --imprecise IMPRECISE

                        Binary value for imprecise (*) hit (can be 1, 0 or

                        0.5), default is 1

  --unconfident UNCONFIDENT

                        Binary value for questionable (?) hit (can be 1, 0 or

                        0.5), default is 0

  --out_prefix OUT_PREFIX

                        Prefix for output file

 

 

テストラン 

オーサーらが準備したテストデータをダウンロードする。

wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR225/ERR225612/ERR225612_*.fastq.gz

Githubにテストデータのリファレンス (S_suis_P17.gbk) とIS配列 (ISSsu3.fasta) もアップされている。それもダウンロードする。

 

中身を確認する。

$ cat ISSsu3.fasta

>ISSsu3

TAGTTAAATGAAACAAAAACAGTACATTTATGATATAATGTATTTATGGCATATTCATTA

GATTTTCGTAAAAAAGTTCTCGCATACTGTGAGAAAACCGGCCGTATTACTGAAGCATCA

GCTATTTTCCAAGTTTCACGTAACACTATCTATCAATGGCTAAAATTAAAAGAGAAAACC

GGCGAGCTTCATCACCAAGTTAAAGGAACCAAGCCAAGAAAAGTGGATAGAGATAAATTA

AAGAATTATCTTGAAACTCATCCAGATGCTTATTTGACTGAAATAGCTTCTGAATTTGAC

TGTCATCCAACAGCTATTCATTACGCCCTCAAAGCTATGGGATATACTCGAAAAAAAAGA

GCTGTACCTACTATGAACAAGACCCTGAAAAAGTAGAACTGTTCCTTAAAGAATTGAATA

ACTTAAGCCACTTGACTCCTGTTTATATTGACGAGACAGGGTTTGAGACATGTTTTCATC

GAGAATATGGTCGCTCTTTGAAAGGTCAGTTGATAAAAGGTAAGGTCTCTGGAAGAAGAT

ACCAGCGGATATCTTTAGTGGCAGGTCTCATAAATGGTGCGCTTATAGCCCCGATGACAT

ACAAAGATACTATGACGAGTGGCTTTTTCGAAGCTTGGTTCAAAATATTCTTACTACCCA

CTTTAGGTAAACCATCTGTTATCATCATGGACAATGCAAAGTTTCATAGGATGAGTAAGC

TAAAAGATTTATGCGAGGAGCAGGGACATAGACTTTTACCACTTCCTCCTTACTCACCGG

AATATAATCCCATTGAGAAAATATGGGCTCACATCAAAAAACACCTCAGAAGAGTATTGC

CAAATTGCGATACTTTTCTTGAGGCACTTTCGTCCTGCTCTTGTTTCAGTTGACTA

 

$ head -40 S_suis_P17.gbk

LOCUS       AM946016             2007491 bp    DNA     circular BCT 14-JUL-2009

DEFINITION  Streptococcus suis P1/7 complete genome.

ACCESSION   AM946016

VERSION     AM946016.1  GI:251819067

DBLINK      BioProject: PRJNA352

KEYWORDS    complete genome.

SOURCE      Streptococcus suis P1/7

  ORGANISM  Streptococcus suis P1/7

            Bacteria; Firmicutes; Bacilli; Lactobacillales; Streptococcaceae;

            Streptococcus.

REFERENCE   1

  AUTHORS   Holden,M.T., Hauser,H., Sanders,M., Ngo,T.H., Cherevach,I.,

            Cronin,A., Goodhead,I., Mungall,K., Quail,M.A., Price,C.,

            Rabbinowitsch,E., Sharp,S., Croucher,N.J., Chieu,T.B., Mai,N.T.,

            Diep,T.S., Chinh,N.T., Kehoe,M., Leigh,J.A., Ward,P.N.,

            Dowson,C.G., Whatmore,A.M., Chanter,N., Iversen,P., Gottschalk,M.,

            Slater,J.D., Smith,H.E., Spratt,B.G., Xu,J., Ye,C., Bentley,S.,

            Barrell,B.G., Schultsz,C., Maskell,D.J. and Parkhill,J.

  TITLE     Rapid evolution of virulence and drug resistance in the emerging

            zoonotic pathogen Streptococcus suis

  JOURNAL   PLoS ONE 4 (7), E6072 (2009)

   PUBMED   19603075

  REMARK    Publication Status: Online-Only

REFERENCE   2  (bases 1 to 2007491)

  AUTHORS   Holden,M.T.G.

  TITLE     Direct Submission

  JOURNAL   Submitted (10-MAR-2008) Holden M.T.G., Pathogen Genomics, Sanger

            Institute Wellcome Trust, Wellcome Trust Genome Campus, Hinxton,

            Cambridge, CB10 1SA, UNITED KINGDOM

FEATURES             Location/Qualifiers

     source          1..2007491

                     /organism="Streptococcus suis P1/7"

                     /mol_type="genomic DNA"

                     /strain="P1/7"

                     /db_xref="taxon:218494"

     gene            1..1374

                     /gene="dnaA"

                     /locus_tag="SSU0001"

                     /gene_synonym="dnaH"

     CDS             1..1374

 

準備ができたらランする。gzip圧縮にも対応している。".fastq.gz”になっている必要がある。

ismap --reads ERR225612_*.fastq.gz --queries ISSsu3.fasta --reference S_suis_P17.gbk --output_dir outdir
  • --reads    Paired end reads for analysing (can be gzipped)
  • --queries   Multifasta file for query gene(s) (eg: insertion sequence) that will be mapped to.
  •   --reference    Reference genome for typing against in genbank format
  •   --output_dir   Location for all output files (default is current directory).
  •   --log    Prefix for log file. If not supplied, prefix will be current date and time.

 

 

テストデータの出力

f:id:kazumaxneo:20170703113539j:plain

gbkファイルを入力すると、このようにIS挿入位置の遺伝子をコールしてくれる。

 

  • ISのクエリ配列は--readsを何度も書くことで複数入力できる。
  • リファレンスはfastaかマルチfasta、またはgenbank形式に対応している。テストデータではgenbank形式のリファレンスを用意している。
  • paired-endシーケンスデータはfastqの非圧縮、またはgzip圧縮に対応している。ただし名前はテストのようにペアとわかる構造の名前でなくてはならない。
  • アセンブルしたcontigや、それにアノテーションをかけgbkにしたファイルを使ってもよい。

 

本手法はゲノムにない新規の挿入でも、配列がわかれば挿入位置を拾ってくれる。バクテリアのゲノムだとランニングタイムも数分で快適に使える。トランスタギングしたような株からIS挿入位置を探したいならまずこのツールを試し、それからペアードエンドデータから手動で探す方法と比較してみるとよいと思われる。

 

* 他の手法と同様に、日本語のフォルダパスの下の方にあると動作しないので気をつける。一度エラーになりました。samtoolsは度々仕様が変わりますが、1.4.1でも動作します。

 

引用

ISMapper: identifying transposase insertion sites in bacterial genomes from short read sequence data

Elizabeth Hénaff, Luís Zapata, Josep M. CasacubertaEmail author and Stephan Ossowski

BMC Genomics. 2015 Sep 3;16:667.