macでインフォマティクス

macでインフォマティクス

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

illuminaと454の前処理ツール seqyclean

 

 最新のハイスループットシーケンス機器は大量のデータを生するが、これにはシーケンスエラー、シーケンスアダプタ、汚染されたリードなどのノイズが含まれていることがよくある。このノイズはゲノミクス研究を複雑にする。シーケンスノイズを減らすために多くの前処理ソフトウェアツールが開発されているが、それらの多くは複数の技術からのデータを扱うことができず、1種類以上のノイズに対処するものはほとんどない。 ここではSeqyCleanという包括的な前処理ソフトウェアパイプラインを紹介する。 SeqyCleanは、ハイスループットシーケンスデータ内の複数のノイズ源を効果的に除去し、著者らのテストによれば、他の利用可能な前処理ツールよりも優れている。 SeqyCleanによる前処理データが最初にde novoゲノム構築とゲノムマッピングの両方を改善することを示す。私たち(著者ら)はアイダホ大学のバイオインフォマティクス・進化論研究所(IBEST)のゲノミクスコアでSeqyCleanを広く使用し、テストデータとproduction dataの両方で検証してきた。 SeqyCleanは、MITライセンス下でhttp://github.com/ibest/seqycleanからオープンソースソフトウェアとして入手できる。

 

SeqyClean offers: (Githubより)

  1. Adapter/key/primers filtering
  2. Vector and contaminants filtering.
  3. Quality trimming.
  4. Poly A/T trimming.
  5. Overlapping paired reads.

 

 

インストール

macos10.14の miniconda3-4.0.5環境でcondaを使って導入した。

依存

Clone or download the repository. Then cd to seqyclean home folder, and type make.

  • zlib
  • make

本体 GIthub

#Bioconda(link)
conda install -c bioconda -y seqyclean

seqyclean -h

$ seqyclean -h

Version: 1.10.09 (2018-10-16)

**********************************************************************************************************************

usage: ./seqyclean libflag input_file_name_1 [libflag input_file_name_2] -o output_prefix [options]

 

Common arguments for all library types:

   -h, --help - Show this help and exit.

   -v <filename> - Turns on vector trimming, default=off. <filename> - is a path to a FASTA-file containing vector genomes.

   -c <filename> - Turns on contaminants screening, default=off, <filename> - is a path to a FASTA-file containing contaminant genomes.

   -k <value> - Common size of k-mer, default=15

   -d - Distance between consecutive k-mers, default=1

   -kc <value> - Size of k-mer used in sampling contaminat genome, default=15

   -qual <max_average_error> <max_error_at_ends> - Turns on quality trimming, default=off. Error boundaries: max_average_error (default=0.01), max_error_at_ends (default=0.01)

   -bracket <window_size> <max_avg_error> - Bracket window_size (default=0.794) and maximum_average_error (default=0.794) for quality trimming

   -window window_size max_avg_error [window_size max_avg_error ...] - Parameters for window trimming. There are two windows with size of 50 and 10bp and max_avg_err of 0.794 by default.

   -ow - Overwrite existing results, default=off

   -minlen <value> - Minimum length of read to accept, default=50 bp.

   -polyat [cdna] [cerr] [crng] - Turns on poly A/T trimming, default=off. Parameters: cdna (default=10) - maximum size of a poly tail, cerr (default=3) - maximum number of G/C nucleotides within a tail, cnrg (default=50) - range to look for a tail within a read.

   -verbose - Verbose output, default=off.

   -detrep - Generate detailed report for each read, default=off.

   -dup [-startdw 10][-sizedw 35] [-maxdup 3] - Turns on screening duplicated sequences, default=off. Here: -startdw (defalt=10) and -sizedw (default=25) are starting position and size of the window within a read, -maxdup (default=3) - maximum number of duplicated sequences allowed.

   -no_adapter_trim - Turns off trimming of adapters, default=off.

Roche 454 only arguments:

   -t <value> - Number of threads (not yet applicable to Illumina mode), default=4.

   -fastq - Output in FASTQ format, default=off.

   -fasta_out - Output in FASTA format, default=off.

   -m <filename> - Using custom barcodes, default=off. <filename> - a path to a FASTA-file with custom barcodes.

Illumina paired- and single-end arguments:

   -1 <filename1> -2 <filename2> - Paired-end mode (see examples below)

   -U <filename> - Single-end mode

   -shuffle - Store non-paired Illumina reads in shuffled file, default=off.

   -i64 - Turns on 64-quality base, default = off.

   -adp <filename> - Turns on using custom adapters, default=off. <filename> - FASTA file with adapters

   -at <value> - This option sets the similarity threshold for adapter trimming by overlap (only in paired-end mode). By default its value is set to 0.75.

   -overlap <value> - This option turns on merging overlapping paired-end reads and <value> is the minimum overlap length. By default the minimum overlap length is 16 base pairs.

   -new2old - Switch to fix read IDs, default=off ( As is detailed in: http://contig.wordpress.com/2011/09/01/newbler-input-iii-a-quick-fix-for-the-new-illumina-fastq-header/#more-342 ).

   -gz - compressed output (GZip format, .gz).

   -alen - Maximum adapter length, default=30 bp.(only for paired-end mode)

Examples

Roche 454:

./seqyclean -454 test_data/in_001.sff -o test/Test454 -v test_data/vectors.fasta

Paired-end Illumina library:

./seqyclean -1 test_data/R1.fastq.gz -2 test_data/R2.fastq.gz -o test/Test_Illumina

Single-end Illumina library:

./seqyclean -U test_data/R1.fastq.gz -o test/Test_Illumina

Please ask Ilya by email: ilya.zhbannikov@duke.edu in case of any questions.

 

 

実行方法

アダプタートリミングを実行する(*1)。

#illumina paired-end(single-endは"-U"を使う)
seqyclean -1 pair_1.fq.gz -2 pair_2.fq.gz -o illumina_output

#454の.sff
seqyclean 454 input.sff -o 454_output -t 8
  • -t <value>    Number of threads (not yet applicable to Illumina mode), default=4.

ペアエンドfastqの出力は同期される。同期機能を無効化するなら"-shuffle"をつける。

 

 

アダプタートリミング、クオリティトリミング、ベクター配列(vector.fasta)のリードフィルタリング、汚染配列(contamination.fasta)のフィルタリングを実行する。

#illumina paired-end
seqyclean -1 pair_1.fq.gz -2 pair_2.fq.gz -o illumina_output \
-qual -v vector.fasta -c contamination.fasta
  • -qual     Turns on quality trimming, default=off
  • -v <filename>     Turns on vector trimming, default=off. <filename> - is a path to a FASTA-file containing vector genomes.
  • -c <filename>     Turns on contaminants screening, default=off, <filename> - is a path to a FASTA-file containing contaminant genomes.

 

 

引用

SeqyClean: A Pipeline for High-throughput Sequence Data Preprocessing

Ilya Y. Zhbannikov, Samuel S. Hunter, James A. Foster, Matthew L. Settles

Conference Paper

ACM-BCB '17 Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology,and Health Informatics
Pages 407-416

 

*1

アダプターP5 / P7 +インデックスI5 / I7 +リンカー、のアダプター配列が除去される。

 

dockerイメージ


アセンブリグラフからプラスミドを検出する HyAsP

 

 プラスミドはバクテリアで一般的なextra-chromosomalのDNA分子である。プラスミドは、それらの長さ(それらはchromosomeよりはるかに短い傾向がある)、コピー数(プラスミドは細胞内に複数のコピーで存在する場合がある)およびGC含有量などの様々な特徴においてchromosomeとは異なる。それらは水平遺伝子導入において、そしてそれ故、ビルレンス因子および抗生物質耐性の伝達において重要な役割を果たす(Carattoli、2013; Dolejska and Papagiannitsis、2018)。それ故、細菌サンプル中のプラスミドの同定は、薬剤耐性細菌の増殖に対する緩和戦略にとって重要である。

 全ゲノムシーケンシング( WGS)データを使用する方法に焦点を合わせて、最近プラスミド検出を行う種々のアプローチが探求されている。最近のプラスミド検出とアセンブリのレビューは (Arredondo-Alonso et al., 2017; Laczny et al., 2017; Orlek et al., 2017)を参照する。既存のツールは、de-novoメソッドとリファレンスベースのメソッドという2種類のアプローチに区別できる。 Recycler(Rozov et al、2017)およびPlasmidSPAdes(Antipov et al、2016)は前者の方法群に属する。Recyclerは、リードデプスおよび長さの特徴に基づいてアセンブリグラフのサイクルを繰り返し剥がすことによってプラスミドを予測する。PlasmidSPAdesは、プラスミドのリードデプスがchromosomeと異なると仮定する。それはchromosomeのリードデプスを推定し、可能性のあるchromosomeのコンティグを取り除き、そして結果として得られる縮小アセンブリグラフの連結成分からプラスミドを予測する。 MOB-recon(Robertson and Nash、2018)は最近のリファレンスベースの方法である。それはリファレンスプラスミドのデータベースならびに既知のレプリコンおよびrelaxasesのコレクションを使用する。コンティグは、リファレンスデータベースに対してマッピングされ、そして推定のプラスミド単位にグループ分けされ、それらはレプリコンまたはrelaxasesなしでそれらの単位を捨てることによってさらに精製される。

 このAPPLICATIONS NOTEでは、WGSアセンブリ内のプラスミドコンティグを識別し、それらを個々のプラスミドにまとめるためのハイブリッドツールで、リファレンスベースの方法とde-novoの方法を組み合わせた新しいツール、HyAsP(プラスミド用ハイブリッドアセンブラ)を紹介する。既知のプラスミド遺伝子とリードデプスやGC含有量などの特性を考慮して、 66のゲノムを含むデータセットについて、HyAsPの予測品質をPlasmidSPAdesおよびMOB-reconと比較し、HyAsPが一般に両方のツールよりも優れていることを示す。
 HyAsPは、既知のプラスミド遺伝子からの情報、リードデプス、およびGC含有量を使用して、アセンブリグラフからプラスミドコンティグを検出、ビニング、およびアセンブリするためのgreedy algorithmである。アルゴリズムアルゴリズム1の擬似コードのワークフローを含む)および実験についての概要を論文の以下に示す。完全な技術詳細は補足資料に記載されている。(以下略)

 

 

インストール

ubuntu16.04にて、condaの仮想環境を作ってテストした(conda create -n hyasp python=3.5.4)。

依存

HyAsP was developed and tested with the following software dependencies:

  • Python (python, version 3.5.4; packages: Bio, math, numpy, os, pandas, random, subprocess, sys)
  • BLAST+ (makeblastdb tblastn and blastn; version 2.6.0)

FastQC (fastqc, version 0.11.5)

  • sickle (sickle, version 1.33)
  • cutadapt (cutadapt, version 1.16)
  • Trim Galore (trim_galore, version 0.4.5_dev)
  • Unicycler (unicycler, version 0.4.5)

   SPAdes (spades.py, version 3.12.0)
   Racon (racon, version 1.3.0)
   Pilon (pilon-1.22.jar, version 1.22)
   SAMtools (samtools, version 1.5)
   Bowtie 2 (bowtie2 and bowtie2-build; version 2.3.3.1)
   Java (java, version 1.8.0_121)

#blast bioconda (link)
conda install -c bioconda -y blast==2.6.0

#fastqc bioconda (link)
conda install -c bioconda -y fastqc==0.11.5

#sickle bioconada (link)
conda install -c bioconda -y sickle

#cutadapt bioconada (link)
conda install -c bioconda -y cutadapt==1.16

#trim galore bioconada (link)
conda install -c bioconda -y trim-galore==0.4.5

#unicycler bioconada (link)
conda install -c bioconda -y unicycler

本体 GIthub

git clone https://github.com/cchauve/hyasp.git
cd hyasp
python setup.py sdist
pip install dist/HyAsP-1.0.0.tar.gz

> hyasp -h

$ hyasp -h

usage: hyasp [-h] {find,create,map,filter} ...

 

positional arguments:

  {find,create,map,filter}

                        task to be performed

    find                find plasmids in the assembly graph

    create              create gene database from collection of plasmids

    map                 map genes to contigs

    filter              filter gene-contig mapping

 

optional arguments:

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

hyasp map -h

$ hyasp map -h

usage: hyasp map [-h] (-f FROM_FASTA | -g FROM_GFA) [-c] [-v]

                 [--makeblastdb MAKEBLASTDB] [--blastn BLASTN]

                 genes_file mapping_file

 

positional arguments:

  genes_file            FASTA file of plasmid genes to be mapped

  mapping_file          output gene-contig mapping

 

optional arguments:

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

  -f FROM_FASTA, --from_fasta FROM_FASTA

                        file containing the contigs (in FASTA format) to which

                        the genes should be matched

  -g FROM_GFA, --from_gfa FROM_GFA

                        file containing the contigs (as part of assembly graph

                        in GFA format) to which the genes should be matched

  -c, --clean           remove temporary files when mapping is done

  -v, --verbose         print more information

  --makeblastdb MAKEBLASTDB

                        path to makeblastdb executable

  --blastn BLASTN       path to blastn executable

> hyasp filter -h

$ hyasp filter -h

usage: hyasp filter [-h] [-i IDENTITY_THRESHOLD] [-l LENGTH_THRESHOLD] [-f]

                    [-v]

                    genes_file mapping filtered_mapping

 

positional arguments:

  genes_file            FASTA file of plasmid genes used to create the mapping

  mapping               gene-contig mapping (BLAST output (default format 6)

                        to be filtered

  filtered_mapping      output file containing the filtered gene-contig

                        mapping

 

optional arguments:

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

  -i IDENTITY_THRESHOLD, --identity_threshold IDENTITY_THRESHOLD

                        minimum identity in a hit to keep it

  -l LENGTH_THRESHOLD, --length_threshold LENGTH_THRESHOLD

                        minimum fraction of query that has to be matched to

                        keep a hit

  -f, --find_fragmented

                        if set, search for fragmented hits, i.e. several short

                        high-quality hits that together satisfy the length

                        threshold

  -v, --verbose         print more information

 

> hyasp find -h

$ hyasp find -h

usage: hyasp find [-h] [-g MIN_GENE_DENSITY] [-k MIN_SEED_GENE_DENSITY]

                  [-l MIN_LENGTH] [-L MAX_LENGTH] [-r MIN_READ_DEPTH]

                  [-d MIN_PLASMID_READ_DEPTH] [-G MAX_GC_DIFF]

                  [-c MAX_INTERMEDIATE_CONTIGS] [-n MAX_INTERMEDIATE_NT]

                  [-s MAX_SCORE] [-w SCORE_WEIGHTS] [-q] [-o OVERLAP_ENDS]

                  [-b BINNING] [-f FANOUT | -p] [-u] [-N] [-v]

                  assembly_graph genes_file gene_contig_mapping output_dir

 

positional arguments:

  assembly_graph        assembly graph (GFA format)

  genes_file            genes database (FASTA format)

  gene_contig_mapping   gene-contig mapping (BLAST output, default format 6)

  output_dir            output directory

 

optional arguments:

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

  -g MIN_GENE_DENSITY, --min_gene_density MIN_GENE_DENSITY

                        minimum gene density of putative plasmids (plasmids

                        with lower gene density are marked as questionable)

  -k MIN_SEED_GENE_DENSITY, --min_seed_gene_density MIN_SEED_GENE_DENSITY

                        minimum gene density of a contig to be considered as a

                        seed

  -l MIN_LENGTH, --min_length MIN_LENGTH

                        minimum length of putative plasmids (shorter plasmids

                        are marked as questionable)

  -L MAX_LENGTH, --max_length MAX_LENGTH

                        maximum length of putative plasmids (longer plasmids

                        are marked as questionable)

  -r MIN_READ_DEPTH, --min_read_depth MIN_READ_DEPTH

                        minimum read depth of contigs to be used in plasmids

  -d MIN_PLASMID_READ_DEPTH, --min_plasmid_read_depth MIN_PLASMID_READ_DEPTH

                        minimum read depth of putative plasmids (plasmids with

                        lower read depth are marked as questionable)

  -G MAX_GC_DIFF, --max_gc_diff MAX_GC_DIFF

                        maximum difference in GC content between contig to be

                        added and current plasmid

  -c MAX_INTERMEDIATE_CONTIGS, --max_intermediate_contigs MAX_INTERMEDIATE_CONTIGS

                        maximum number of contigs between gene-containing

                        contigs in a plasmid

  -n MAX_INTERMEDIATE_NT, --max_intermediate_nt MAX_INTERMEDIATE_NT

                        maximum sum of lengths of contigs between gene-

                        containing contigs in a plasmid

  -s MAX_SCORE, --max_score MAX_SCORE

                        maximum score for eligible extensions

  -w SCORE_WEIGHTS, --score_weights SCORE_WEIGHTS

                        weights of the score components (comma-separated list

                        of elements <component>=<weight>)

  -q, --keep_subplasmids

                        do not mark plasmids contained by others as

                        questionable

  -o OVERLAP_ENDS, --overlap_ends OVERLAP_ENDS

                        minimum overlap between ends of plasmid sequence for

                        circularisation

  -b BINNING, --binning BINNING

                        number standard deviations the read depth and GC

                        content of plasmids are allowed to differ from the

                        centre of their bin (binning is activated by setting

                        the parameter to any value unequal NaN)

  -f FANOUT, --fanout FANOUT

                        maximum number of predecessors / successors of any

                        contig in a plasmid (resp. contig collection)

  -p, --probabilistic   activates probabilistic extensions

  -u, --use_median      determine average read depth of a plasmid using the

                        median instead of the mean

  -N, --use_node_based  activates node-based extensions

  -v, --verbose         print more information

(hyasp) kazu@kazu:~/Document/hyasp/example/example_files$ 

 

>hyasp_pipeline -h

$ hyasp_pipeline -h

usage: hyasp_pipeline [-h] [-1 FIRST_SHORT_READS] [-2 SECOND_SHORT_READS]

                      [-s SINGLE_SHORT_READS] [-l LONG_READS]

                      [-u UNICYCLER_MODE] [-f] [-n]

                      [--read_qual_threshold READ_QUAL_THRESHOLD]

                      [--min_read_length MIN_READ_LENGTH]

                      [--unicycler_verbosity UNICYCLER_VERBOSITY]

                      [--unicycler_keep UNICYCLER_KEEP]

                      [--identity_threshold IDENTITY_THRESHOLD]

                      [--length_threshold LENGTH_THRESHOLD]

                      [--min_gene_density MIN_GENE_DENSITY]

                      [--min_seed_gene_density MIN_SEED_GENE_DENSITY]

                      [--min_length MIN_LENGTH] [--max_length MAX_LENGTH]

                      [--min_read_depth MIN_READ_DEPTH]

                      [--min_plasmid_read_depth MIN_PLASMID_READ_DEPTH]

                      [--max_gc_diff MAX_GC_DIFF]

                      [--max_intermediate_contigs MAX_INTERMEDIATE_CONTIGS]

                      [--max_intermediate_nt MAX_INTERMEDIATE_NT]

                      [--max_score MAX_SCORE] [--score_weights SCORE_WEIGHTS]

                      [--keep_subplasmids] [--overlap_ends OVERLAP_ENDS]

                      [--binning BINNING] [--fanout FANOUT] [--probabilistic]

                      [--use_median] [--use_node_based] [--verbose]

                      [--fastqc FASTQC] [--sickle SICKLE]

                      [--cutadapt CUTADAPT] [--trim_galore TRIM_GALORE]

                      [--unicycler UNICYCLER] [--spades SPADES]

                      [--racon RACON] [--pilon PILON] [--samtools SAMTOOLS]

                      [--bowtie2 BOWTIE2] [--bowtie2_build BOWTIE2_BUILD]

                      [--java JAVA] [--makeblastdb MAKEBLASTDB]

                      [--tblastn TBLASTN] [--blastn BLASTN]

                      analysis_dir genes_file

 

positional arguments:

  analysis_dir          directory into which all data and outputs are stored

  genes_file            (plasmid) genes to map to contigs

 

optional arguments:

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

  -1 FIRST_SHORT_READS, --first_short_reads FIRST_SHORT_READS

                        FASTQ file of first paired reads (if any)

  -2 SECOND_SHORT_READS, --second_short_reads SECOND_SHORT_READS

                        FASTQ file of second paired reads (if any)

  -s SINGLE_SHORT_READS, --single_short_reads SINGLE_SHORT_READS

                        FASTQ file of single / unpaired reads (if any)

  -l LONG_READS, --long_reads LONG_READS

                        FASTQ file of long reads (if any)

  -u UNICYCLER_MODE, --unicycler_mode UNICYCLER_MODE

                        Unicycler assembly mode (conservative, normal or bold)

  -f, --use_fastqc      if set, FASTQ reads are analysed with FastQC before

                        and after preprocessing

  -n, --no_preprocessing

                        if set, the FASTQ reads are not preprocessed before

                        the assembly

  --read_qual_threshold READ_QUAL_THRESHOLD

                        threshold for trimming low-quality ends

  --min_read_length MIN_READ_LENGTH

                        minimum length of reads after quality / adapter

                        trimming

  --unicycler_verbosity UNICYCLER_VERBOSITY

                        Unicycler level of stdout and log file information (0

                        - 3)

  --unicycler_keep UNICYCLER_KEEP

                        Unicycler level of file retention (0 - 3)

  --identity_threshold IDENTITY_THRESHOLD

                        filtering the gene-contig mapping: minimum identity in

                        a hit to keep it

  --length_threshold LENGTH_THRESHOLD

                        filtering the gene-contig mapping: minimum fraction of

                        query that has to be matched to keep a hit

  --min_gene_density MIN_GENE_DENSITY

                        greedy algorithm: minimum gene density of predicted

                        plasmids (are marked as questionable otherwise)

  --min_seed_gene_density MIN_SEED_GENE_DENSITY

                        greedy algorithm: minimum gene density of a contig to

                        be considered as a seed

  --min_length MIN_LENGTH

                        greedy algorithm: minimum length of predicted plasmid

                        (shorter ones are marked as questionable)

  --max_length MAX_LENGTH

                        greedy algorithm: maximum length of predicted plasmid

                        (longer ones are marked as questionable)

  --min_read_depth MIN_READ_DEPTH

                        greedy algorithm: minimum read depth of contig to be

                        used in plasmids

  --min_plasmid_read_depth MIN_PLASMID_READ_DEPTH

                        greedy algorithm: minimum read depth of predicted

                        plasmids (are marked as questionable otherwise)

  --max_gc_diff MAX_GC_DIFF

                        greedy algorithm: maximum difference in GC content

                        between contig to be added and current plasmid

  --max_intermediate_contigs MAX_INTERMEDIATE_CONTIGS

                        greedy algorithm: maximum number of contigs between

                        gene-containing contigs in a plasmid

  --max_intermediate_nt MAX_INTERMEDIATE_NT

                        greedy algorithm: maximum sum of lengths of contigs

                        between gene-containing contigs in a plasmid

  --max_score MAX_SCORE

                        greedy algorithm: maximum score for potential

                        extensions

  --score_weights SCORE_WEIGHTS

                        greedy algorithm: weights of the score components

                        (comma-separated list of elements

                        <component>=<weight>)

  --keep_subplasmids    greedy algorithm: mark plasmids contained by others as

                        questionable

  --overlap_ends OVERLAP_ENDS

                        greedy algorithm: minimum overlap between ends of

                        plasmid sequence for circularisation

  --binning BINNING     greedy algorithm: setting a threshold value activates

                        plasmid binning

  --fanout FANOUT       greedy algorithm: maximum degree of branching of

                        contigs (per end)

  --probabilistic       greedy algorithm: choose extension probabilistically

  --use_median          greedy algorithm: determine average read depth of a

                        plasmid using the median instead of the mean

  --use_node_based      greedy algorithm: activates node-based extensions

  --verbose             greedy algorithm: prints detailed information during

                        the extension process

  --fastqc FASTQC       path to FastQC executable

  --sickle SICKLE       path to sickle executable

  --cutadapt CUTADAPT   path to cutadapt executable

  --trim_galore TRIM_GALORE

                        path to Trim Galore executable

  --unicycler UNICYCLER

                        path to Unicycler executable

  --spades SPADES       path to SPAdes executable

  --racon RACON         path to Racon executable

  --pilon PILON         path to Pilon executable

  --samtools SAMTOOLS   path to SAMtools executable

  --bowtie2 BOWTIE2     path to Bowtie 2 executable

  --bowtie2_build BOWTIE2_BUILD

                        path to bowtie2_build executable

  --java JAVA           path to Java executable

  --makeblastdb MAKEBLASTDB

                        path to makeblastdb executable

  --tblastn TBLASTN     path to tblastn executable

  --blastn BLASTN       path to blastn executable

 

データベースの準備

cd databases 
hyasp create ncbi_database_genes.fasta -p plasmids.csv \
-b ncbi_blacklist.txt -d -l 500 -m 100 -t GenBank \
-r 2015-12-19T00:00:00Z

ncbi_database_genes.fastaができる。これをラン時に指定する(以下はgenes.fastaと記載)。

 

実行方法

1、アセンブリグラフ(GFA)からplasmidを探索。データベースncbi_database_genes.fastaアセンブリグラフのGFAを指定する。 

hyasp map genes.fasta gcm.csv -g assembly.gfa

結果はblastテキスト出力ライクな形式になる。

 

2、必要に応じて出力をフィルタリングする。

hyasp filter genes.fasta gcm.csv filtered_gcm.csv \
--identity_threshold 0.95 \
--length_threshold 0.95
  • --identity_threshold   identity of hits retained in the mapping (default: 0.95).
  • --length_threshold     Minimum fraction of query (gene) that has be matched to keep a hit (default: 0.95).

 

3、プラスミド配列検出。1の出力か2のフィルタリングした出力を指定する。

hyasp find assembly.gfa genes.fasta gcm.csv output_dir

上に載せたヘルプにあるように、感度調整のパラメータが多いので注意する。

出力

> ls -alth

-rw-rw-r-- 1 kazu kazu 4.2M  7月 13 22:44 tagged_assembly.gfa

drwxrwxr-x 2 kazu kazu 4.0K  7月 13 22:44 .

-rw-rw-r-- 1 kazu kazu  50K  7月 13 22:44 questionable_plasmid_contigs.fasta

-rw-rw-r-- 1 kazu kazu  89K  7月 13 22:44 putative_plasmid_contigs.fasta

-rw-rw-r-- 1 kazu kazu  90K  7月 13 22:44 putative_plasmids.fasta

-rw-rw-r-- 1 kazu kazu  50K  7月 13 22:44 questionable_plasmids.fasta

-rw-rw-r-- 1 kazu kazu  270  7月 13 22:44 contig_chains.csv

推定piasmid配列FASTAの他に、元のGFAファイルに対して推定プラスミドコンティグの色およびラベル情報が追加されたGFAファイルが出力される。

f:id:kazumaxneo:20190713225617p:plain

少なくとも1つの推定プラスミド中に存在するコンティグは青色がアサインされ、1つまたは複数の疑わしいプラスミド中に存在すると水色がアサインされる。

 

 

NGSのリードを指定して、rawリードから直接実行することもできる。

hyasp_pipeline output_dir ncbi_database_genes.fasta \
-1 pair_1.fq -2 pair_2.fq
  • -1    FASTQ file of first paired reads (if any)
  • -2    FASTQ file of second paired reads (if any)
  • -s     FASTQ file of single / unpaired reads (if any)
  • -l      FASTQ file of long reads (if any)

 

シングルエンドショートリードとロングリードを指定。

hyasp_pipeline output_dir ncbi_database_genes.fasta \
-s short_reads.fastq -l long_reads.fastq

 

テスト時はunicyclerを使ったアセンブリのステップでエラーになって結果は確認できなかった。

引用

HyAsP, a greedy tool for plasmids identification
Robert Müller, Cedric Chauve
Bioinformatics, Published: 22 May 2019

 

 

 

 

 

シングルの配列やメタゲノムのbinned.fastaのtaxonomic classificationを行う BASTA

2019 7/13 説明修正

2019 8/1 説明追記

2020 1/21 インストール手順修正

2020 2/4 データベースダウンロード手順修正

2020 4/17 コマンド修正

2020 4/19 binned fastaを使う手順追記

 

 DNAシーケンシング、例えばアンプリコン、メタゲノムおよび全ゲノムシーケンシングは、微生物学および生態学から医学まで、ライフサイエンスの多くの分野において標準的な手順となっている。これらのシークエンシングプロジェクトの大部分における重要なステップは、例えば全ゲノムシークエンシング(WGS)プロジェクトにおける培養汚染を同定するため、またはメタゲノムアセンブリからの特定の連続した配列(コンティグ)を特定の分類学に結び付けるための遺伝物質のtaxonomy分類ステップである。

 長年にわたり、アンプリコンシークエンシングとメタゲノムにおけるNGSリードの解析(Caporaso et al、2010; Schloss et al、2009)に特に焦点を当てて、シークエンシングデータをtaxonomyに分類するためのさまざまなアプローチが開発されている(Darling et al、2014; Huson、Auch、Qi、&Schuster、2007; Wood&Salzberg、2014)。ただし、メタゲノムアセンブラとビニングアルゴリズムの精度の向上、およびOxford NanoporeやPacBioシーケンスなどの第3世代シーケンシング(3GS)テクノロジのリード長の向上には、分類学的な分類のためのより一般的なアプローチが必要になる。 

 Last Common Ancestor (LCA)のtaxonomy推定は、最小共通祖先(lowest common ancestor)(Wood&Salzberg、2014)またはコンセンサス(Caporaso et al、2010)分類法としても知られており、NGSリードのtaxonomy分類に広く使用されているアプローチである。 (Wood&Salzberg、2014)は、長さkの配列を共有するゲノム、いわゆるk-merを、クエリを用いて同定することによって、NGSメタゲノムデータから微生物群集プロファイルを推定した。読み取られたクエリのtaxonomyは、同じk-merを含むすべてのゲノムのLCAとして推定される。同様に、メタゲノム解析ツールキットMEGAN(Huson et al。、2007)は、配列類似性検索ツールDIAMOND(Buchfink、Xie、&Huson、2015)に基づくLCA推定を使用して、NGSリードの分類と機能プロファイルをリンクさせている。 LCA推定の別の例は、アンプリコン配列の分類のためにリードの3つの最良のヒットを使用する、アンプリコン分析フレームワークQIIME(Caporaso et al., 2010)において実施されたBLASTベースのtaxonomy推定である。 LCAアルゴリズムは、十分に精緻化されたターゲットデータベースに依存しているが、非常に正確であることが示されており(McIntyre et al、2017)、使用する比較ツールによっては計算効率が高くなる。しかし、現在のLCA実装はNGSリードまたは短い配列に限定されており、予測して得たアミノ酸配列、ゲノムおよびメタゲノムプロジェクトからアセンブリされたコンティグ、または3GS技術によって生み出されるロングリード配列などを分類する能力がない。さらに、著者らの知る限りでは、既存のLCA実装は、事前定義されたワークフローの外側でシーケンス分類を提供しない、より大きなソフトウェアフレームワークに組み込まれている。

 ここでは、塩基配列分類アノテーションツールBASTAを提示する。これは、長さやシーケンシングプラットフォームに関係なく、LCAアルゴリズムヌクレオチドおよびアミノ酸配列に拡張する高度にカスタマイズ可能なLCA実装を提供する。 BASTAは、多くの一般的な配列類似性検索ツールと一緒に使用することができる、e.g. BLAST (Altschul, Gish, Miller, Myers, & Lipman, 1990)、Diamond (Buchfink et al., 2015)。そしてBASTAは単一のシーケンスならびにシーケンスグループおよびシーケンスビンにアノテーションを付けることができる。 BASTAはpython 2.7で書かれ、GNU General Public Licenseの下で開発されたコマンドラインツールである。すなわちオープンソースhttps://github.com/timkahlke/BASTAで無料で入手できる。インストールを簡単にするために、BASTAはパッケージマネージャCondaを使ってインストールすることも、pythonのsetuptoolsを使ってネイティブインストールすることもできる。詳細なインストール手順とドキュメントはhttps://github.com/timkahlke/BASTA/wikiにある。
 BASTAは、sequence_id - > taxon_id(ftp://ftp.ncbi.nih.gov/pub/taxonomy)の形式の対応するマッピングファイルと共に、NCBIによって公表されたゲノムのtaxonomyを利用している。ヒットしたtaxonomiesは、ヒットのアクセッション番号をそのNCBI tax IDにマッピングすることによって決定される。その後、それは7レベルのtaxonomyとしてresolveされる(論文図1)。 BASTAをインストールしたら、コアtaxonomyをダウンロードしてlevelDBデータベースにインポートする必要がある。このプロセスは完全に自動化されており、BASTAのtaxonomyコマンドに基づいて実行される。 BASTAのダウンロードコマンドは、RefSeq(O'Leary et al、2016)、Uniprot(The Uniprot Consortium、2017)、protein data bank(pdb)( Berman et al、2000)など、最も一般的なシーケンスデータベースのローカルデータベースへのマッピングの自動ダウンロードとそれに続くインポートを実装している。さらに、BASTAのcreate_dbコマンドを使用してカスタムマッピングファイルをインポートすることもできる。(以下略)

 

HP

https://uts-c3.github.io/project/basta/

wiki

https://github.com/timkahlke/BASTA/wiki

 

f:id:kazumaxneo:20190708202756p:plain

BASTA workflow. 論文より転載

 

インストール

ubuntu16.0.4のminiconda3.4.0.5環境でテストした。

依存

ASTA requires the non-standrad python packages.

本体 Github

#ここではcondaの仮想環境を作って導入, diamondも使うので入れる。
conda create -n basta -y
conda activate basta
conda install -c bioconda -c bnoon -c timkahlke basta diamond

テスト

python -m unittest discover basta

# python -m unittest discover basta

.........--2019-07-08 20:31:48--  http://www.google.com/index.html

Resolving www.google.com... 172.217.26.100, 2404:6800:400a:808::2004

Connecting to www.google.com|172.217.26.100|:80... connected.

HTTP request sent, awaiting response... 200 OK

Length: unspecified [text/html]

Saving to: '/root/.pyenv/versions/miniconda3-4.0.5/lib/python2.7/site-packages/basta/index.html'

 

/root/.pyenv/versions/miniconda3-4.0.5/lib     [ <=>                                                                                   ]  12.33K  --.-KB/s    in 0.001s 

 

2019-07-08 20:31:48 (12.7 MB/s) - '/root/.pyenv/versions/miniconda3-4.0.5/lib/python2.7/site-packages/basta/index.html' saved [12623]

 

........

----------------------------------------------------------------------

Ran 17 tests in 0.127s

 

OK

 

O.k ! 

basta --help

# basta --help

usage: basta [-h] {sequence,single,multiple,download,create_db,taxonomy} ...

 

Basic sequence taxonomy assignment

 

positional arguments:

  {sequence,single,multiple,download,create_db,taxonomy}

                        sub-command help

 

optional arguments:

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

> basta download -h

$ basta download -h

usage: basta download [-h] [-d DIRECTORY] [-f FTP]

                      {wgs,prot,est,gss,gb,pdb,uni}

 

Download NCBI taxonomy file(s)

 

positional arguments:

  {wgs,prot,est,gss,gb,pdb,uni}

                        Type of mapping file to be downloaded (prot, est, wgs,

                        gss or gb)

 

optional arguments:

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

  -d DIRECTORY, --directory DIRECTORY

                        Directory of mapping files (default:

                        $HOME/.basta/taxonomy)

  -f FTP, --ftp FTP     URL to NCBI ftp for accession mapping (default:

                        ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/)

  

データベースの準備 

#NCBI taxonomy databaseのセットアップ(7 level taxonomy file)
basta taxonomy -d download_dir

#NCBI's protein-to-taxonIDとgenbank-to-taxonID mapping fileのセットアップ(link
#下で使うデータベースに合わせたダウンロードする。estを実行するならbasta download est
#ここではprotとgbをダウンロード
basta download gb -d download_dir
basta download prot -d download_dir

以下の既知データベースが利用できる(wikiより)。 

  • prot - protein-to-taxonID mapping file (protein sequences hosted at the NCBI)
  • uni - uniprot-to-taxonID mapping file (complete uniprot, will also be imported into database prot)
  • gb - genbank-to-taxonID mapping file (for most nucleotide databases)
  • wgs - whole-genome-sequence-to-taxonid mapping file
  • pdb - protein-database-to-taxonID mapping file
  • est - est-to-taxonID mapping file
  • gss - gss-to-taxonID mapping file

保存場所を指定しない場合、defaultではBASTA_INSTALL_DIR/taxonomy/に保存される。

 

実行方法

様々な使い方があるが、典型的な使い方の1つは、クエリのプロテイン配列とUniprot databaseのproteomeとのblast結果を使ったTaxonomic classificationである。

 

 1、uniprot databaseに対してblast検索を実行する。Diamondを使う。例えばUniProt KnowledgebaseのUniref90をダウンロードしてdiamond makedbでデータベースにしておく。

Diamondを使いblast検索を実行(Uniref90だとかなり重たいので注意、用途に合わせて変更する)。

diamond blastx --query input.fa \
--db uniprot_ref_proteomes.diamond.dmnd \
--outfmt 6\
--sensitive \
> blast.out

 

2、bastaを実行。<MAPPING_FILE_TYPE> はデータによって変更する。"basta sequence"だと配列それぞれのLCA推定、"basta single"だと全ヒットから1つのLCA推定になる。

#ここでは<MAPPING_FILE_TYPE> は”prot”にする(上のデータベースの準備参照)。
#basta sequence
basta sequence -d download_dir blast.out output prot

#basta single
basta single -d download_dir blast.out output prot
  • MAPPING_FILE_TYPE = one of either prot, gb, est, gss,pdb or wgs 

E.coliのproteomeを使った時のsingleの出力。

$ cat output 

Sequence Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

 

E.coliのproteomeを使った時のsequenceの出力。

> head output

$ head output 

NP_414594.1 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414596.1 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414611.1 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414671.1 Unknown

NP_414695.1 Unknown

NP_414700.1 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414705.1 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414756.2 Bacteria;Proteobacteria;Gammaproteobacteria;Enterobacterales;Enterobacteriaceae;

NP_414805.1 Unknown

NP_414815.1 Unknown

 

 

メタゲノムのbinned fastaなど複数のblast hitファイルからLCA推定を行うなら"basta multiple”を使う。

#!/bin/sh
#prooka annotation
mkdir blast_outdir

for file in `\find *fa -maxdepth 1 -type f`; do
genome=${file}
out=${file%.fa}
diamond blastx --query $genome --db uniprot_ref_proteomes.diamond.dmnd --outfmt 6 --sensitive > blast_outdir/${out}
done

 basta_loop.shとして保存し、binned_fastaディレクトリで実行する。

./ basta_loop.sh

 

 得られたblast_outdirを指定してbasta multipleを実行。

basta multiple blast_outdir/ basta_out gb

 

引用

BASTA – Taxonomic classification of sequences and sequence bins using last common ancestor estimations

Tim Kahlke Peter J. Ralph

 


 

de novo transcriptome解析のクラスタリングとclosely rellatedな種の情報を用いたアノテーションを行う Grouper

 

 シーケンシング技術の進歩により、モデル生物の範囲を超えてトランスクリプトームを効率的かつ正確に探索することが可能になった(Ekblom and Galindo、2011; Marioni et al、2008)。トランスクリプトームシークエンシングは、高品質のリファレンスゲノムを持たない種でも遺伝子発現を理解するための扉を開く。これは既知の種の大多数を構成している(Martin and Wang、2011)。トランスクリプトームシークエンシングのいくつかの一般的な用途には、変異検出、フュージョンおよび選択的スプライシングイベントの発見、ならびにdifferential expression analysisが含まれる(Soumana et al、2015; Stubben et al、2014)。転写物の配列が先験的に知られていない生物のショートリードシーケンスでは、決定的な最初のステップは、ショートシークエンシングリードを使ってトランスクリプトームをアセンブリすることである。
 トランスクリプトームは典型的にはディープにシーケンシングされ、その結果、数千から数億のリードが生じ、それを元の転写産物シーケンスを再構築するために de novo assemblyというプロセスを通してアセンブリする必要がある。これらの転写産物は後の分析のための「リファレンス」として後で機能するであろう。例えば、differential expression分析パイプラインでは、リードが最初にアセンブリされたトランスクリプトマッピングされ、トランスクリプトの量が推測される(Li and Dewey、2011)。
Trinity(Haas et al、2013)、Oasis(Schulz et al、2012)、Trans-ABySS(Robertson et al、2010)など、冗長性​​とリード間の重複を利用してショートリードから完全長のトランスクリプトアセンブリするための多くの一般的なツールがある。同様に、これらの方法で生成されたアセンブリの品質を向上させるために使用できるツールがある(Cabau et al、2017; Durai and Schulz、2016)。ただし、強力なアルゴリズムヒューリスティックが採用されているにもかかわらず、最終的な出力シーケンス(つまりコンティグ)は、完全長のトランスクリプトを表していないことがよくある。言い換えれば、多くの転写産物は重複しないコンティグのセットに断片化されている。これは、シーケンシングされたリードにおけるエラー、オルタナティブスプライシングおよびパラロガス遺伝子から生じるトランスクリプトーム複雑性、シーケンシングされた分子の長さにわたる不均一なまたは不十分なカバレージアセンブリ(例えば、解決される計算上の問題が根底にある生物学とうまく一致しない)ならびに以下に使用される基礎的アプローチの欠陥を含む。これらのすべてが、アッセイされた実際の転写産物のセットよりもはるかに大きい出力コンティグ集団をもたらし得る。
 潜在的に低い品質のde novoアセンブリアウトプットのために、転写産物レベル(すなわち、コンティグレベル)分析の代わりに遺伝子レベルのdifferential expression分析を追求することがより有望でありそして頑強である。遺伝子レベルの情報を得るためには、同じ転写産物の一部を表すコンティグと同じ遺伝子のアイソフォームを表すコンティグを推測し、それらをグループ化する必要がある。これは、Davidson and Oshlack(2014)、Ptitsyn et al(2015)、Srivastava et al(2016)が提案しているものである。 Corsetは、コンティグの共有カウントに基づいて階層的クラスタリングを計算し、RapClustは、コンティグの存在量に基づいて重み付けされた、リードを共有するエッジとコンティグを表すmapping ambiguity graphと呼ばれるスパースグラフをクラスタリングした。さらに、これらのクラスターの意味のある解釈と様々な分析のために、de novoアセンブリのコンティグが何を表しているかについてのいくつかの概念を持つことが重要である(Garber et al、2011)。多くのケースで、我々は研究中の非モデル生物とclosely relatedな種からのゲノムまたはトランスクリプトームを使いアノテーションを付ける。この情報を活用して、de novoアセンブリのコンティグに正確にアノテーションを付け、differential expressionな遺伝子(DEGs)の機能に関する理解を深めることができる。伝統的に、BLASTのバリエーション(Altschul et al、1990)を用いてこのアノテーション付けを行い、次いで遺伝子オントロジー(GO)analysisを完成させる(Ji et al、2012; Parchman et al、2010)。
 現在、下流分析を改善し、そのような分析からの結果の解釈を容易にするために、新たなアセンブリを処理するいくつかの方法が存在する。しかし、それぞれの方法には限界があり、Corset以外には、クラスタリングや新規のコンティグへのアノテーション付けのための完全なパイプラインを提供するツールは存在しない。 RapClustの枠組みを包み込んで改良し、de novoトランスクリプトームアセンブリの正確で効率的な処理に使用できるGrouperを紹介する。Grouperの基礎となるアルゴリズムはde novoアセンブリからグラフに情報を転送するため、アセンブリされたコンティグの直接解析を超え、将来のde novoトランスクリプトーム解析の有用な方向性を示唆している。完全なGrouperパイプラインは、differential expression analysesに使用できる定量情報も提供する。 Grouperソフトウェアはhttps://github.com/COMBINE-lab/grouperにて自由に利用できる。

 

 

  Grouperには、クラスタリングモジュールとラベリングモジュールの2つの主要モジュールがある。 前者はツールRapClustに基づいており、迅速な転写産物レベル定量のためSailfishまたはSalmonツールの下流で実行するように設計されている。 Grouperのラベリングモジュールは、closely relatedな種のゲノムのアノテーション情報を使い、de novoアセンブリ内のコンティグにアノテーションを付けるために使用できる。 

 

インストール

ubuntu16.04にて、condaの仮想環境を作ってテストした(conda create -n Grouper python=2.7.15)。

依存

  • python2
  • Click
  • PyYAML
  • Pandas
  • NumPy
  • Networkx v1.11
  • mcl

本体 Github

#依存も含めpipで導入可能
pip install biogrouper

#MCLがなかったので別途導入した
conda install -c bioconda -y mcl

> Grouper -h

$ Grouper -h    

Usage: Grouper [OPTIONS]

 

Options:

  --config TEXT  Config file describing the experimental setup  [required]

  -h, --help     Show this message and exit.

 dockerイメージ(未テスト)

docker pull combinelab/grouper

 

実行方法

1、Githubの解説の通り、ヒトRNA seqのバイオプロジェクトPRJNA162905 (link) の6つのSRAファイルをダウンロードする。

ここではfasterq-dumpを使っている(紹介)。forで回す。

sample=(SSRR493366 SRR493367 SRR493368 SRR493369 SRR493370 SRR493371) 

for name in ${sample[@]}; do
fasterq-dump $name -e 8 -p
done

 

 

2、リファレンスtranscriptsのindexing(NCBIよりGRCh38のRefSeq Transcriptsダウンロードして使用 => link

salmon index -p 20 -t GRCh38_latest_rna.fna.gz -i ref_index
  • -p  [ --threads ] arg (=2) Number of threads to use (only used for computing bias features)
  • -t  [ --transcripts ] arg Transcript fasta file.
  • -k  [ --kmerLen ] arg (=31) The size of k-mers that should be used for the quasi index.
  • -i  [ --index ] arg Salmon index.

ref_index/が出力される。これをstep3で指定する。 

 

 

3、salmonによる定量.-dumpEqと--writeOrphanLinksオプション付きで実行する。ここではGNU parallelで6ジョブを並列処理している。

parallel -j 6 "samp={}; salmon quant -i ref_index -l a \
-1 reads/{$samp}_1.fastq -2 reads/{$samp}_2.fastq \
-o {$samp}_quant --dumpEq --writeOrphanLinks -p 4" \
::: SRR493366 SRR493367 SRR493368 SRR493369 SRR493370 SRR493371

 > ls -alth

$ ls -alth

合計 111M

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:29 SRR493371_quant

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:29 SRR493368_quant

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:29 SRR493370_quant

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:28 SRR493367_quant

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:28 SRR493369_quant

drwxrwxr-x  5 kazu kazu 4.0K  7月 10 18:28 SRR493366_quant

 

それぞれのディレクトリ({$samp}_quant)に各サンプルの定量結果が出力される。

 

 

4、Grouperによるクラスタリング

 yamlファイルを準備して、そこにグループ、step3で得た処理群とコントロール群それぞれの定量ファイルパス、出力ディレクトリなどの情報を記載する。Githubの通りだと、以下のようになる。

conditions:

    - Control

    - HOXA1 Knockdown

samples:

   Control:

       - SRR493366_quant

       - SRR493367_quant

       - SRR493368_quant

   HOXA1 Knockdown:

       - SRR493369_quant

       - SRR493370_quant

       - SRR493371_quant

outdir: human_grouper

orphan: True

mincut: True

 

ファイル名"config.yaml"として保存し、Grouperを実行。

Grouper --config config.yaml

このステップにはかなりの時間がかかる(*1)。ランが終わると、指定したhuman_grouper/に複数ファイルが出力される。mag.flat.clustがクラスター情報のファイル。 tximportのようなツールと互換性のある "transcript-to-gene"マッピング形式で計算されたクラスター情報が含まれている。

 

近縁種のtranscriptsのアノテーションは、transcripts.fastaを用意し、それをyamlファイルに記載することで機能します。詳細はGithubで確認してください。

 

引用

Grouper: graph-based clustering and annotation for improved de novo transcriptome analysis
Laraib Malik, Fatemeh Almodaresi, Rob Patro
Bioinformatics, Volume 34, Issue 19, 01 October 2018, Pages 3265–3272

 

関連


*1

テスト時は数日かかった。

 

 

Pacbioシーケンシングリードのオーバーラップ検出感度を改善する GroupK

 

 リード長の増加により、第3世代のシークエンシングでゲノムアセンブリのギャップを埋め[ref.1, 2]、構造の変化を明らかにし[ef.13]、トランスクリプトームシークエンシングで遺伝子アイソフォームをより正確に定量できるようになった[ef.14]。さらに、ロングリードを使用することは、微生物群集構造を明らかにし、微生物群集の種内変動を解読することにおいて有望である[ref.5、6]。

 第三世代のシークエンシングデータを用いたゲノムの構築には、専用の方法とツールが必要になる。既存のゲノム構築ツールは、主に2つのタイプのグラフモデル、すなわちオーバーラップグラフおよびde Brujinグラフを利用する。エラー率が低い場合、de Bruijnグラフには、シーケンスサイズが大きくなってもグラフサイズが大きく増加しないという理論的な利点がある。これは通常、イルミナのデータセットでは高い値である。第3世代のシークエンシングデータでは、高いエラー率と低いカバレッジにより、オーバーラップグラフがゲノム構築に適した選択肢となる[ref.7]。オーバーラップグラフを構築する際の重要なステップは、オーバーラップを共有するリードペアを同定することであり、これはそれらリードが基礎となるゲノム内の同じ遺伝子座からシーケンシングされたことを示す。オーバーラップアラインメントを実施するために利用可能な多数のシーケンスアラインメントプログラムがあるが[ref.8、9]、それらの大部分は動的計画法に頼っており、そしてハイスループットシーケンシングデータにとって計算コストがかかりすぎる。エラー率が高いため、BWT(Burrows-Wheeler変換)またはハッシュテーブル[ref.10、11]を使用した既存のショートリードのオーバーラップ検出ソフトウェアは、ロングリードに直接適用することはできない。

 エラーが発生しやすいロングリードの重複を検出するために、現在2つの方法が採用されている。1つの戦略は、オーバーラップ検出の前に、Pacbio(パシフィックバイオサイエンス)およびONT(オックスフォードナノポア)データにおけるシーケンシングエラーを修正することを試みるものである。様々なシーケンシングエラー訂正ツールが存在する[ref.7、12]。それらのうちのいくつかは、少なくとも2つのシーケンシングライブラリーおよび数種類のシーケンシングラン調製を必要とし、したがって多くの用途にとって費用効果が高くないハイブリッドシーケンシングに依存している。他のものはロングリードだけを使用してエラー訂正を行う。代表的な方法の1つが、ChinらのHGAP [ref.12]に記載されており、その性能は、より高いリードカバレッジ範囲で向上する。 HGAPのようなアライメントベースのエラー訂正方法にとって、重要なステップは迅速にアライメントできるリードを識別することである。本質的に、オーバーラップ検出に使用される技術はアライメント検出にも同様に使用することができる。

 2番目のカテゴリは、エラー訂正の難しさを回避し、rawリードを使用して重複を識別する。 PacBioおよびONTデータには、さまざまな近似類似検索方法が適用されている[ref.13]。これらは一般的にシードチェーンアラインメントの手順に従う[ref.14]。シードベースのフィルタリングステップは、感度と計算効率の間のトレードオフを制御する上で重要な役割を果たす。通常、これらの方法では、フィルタリングステップとして短い文字列の(完全)一致を使用する。短い文字列またはk-merの一致には、2つのシーケンス間のk個の連続した文字の完全一致が必要である。直感的には、オーバーラップするリードは、オーバーラップしないリードよりも一般的なk-merを共有する傾向がある。したがって、共有k-mer数を迅速に見つけることができる戦略を適用することができる。このセクションでは、いくつかの最先端のオーバーラップ検出ツールの主な戦略をまとめる。次のセクションで、我々(著者ら)の方法と既存の方法との違いを強調する。

(3段落省略)

 この研究では、ゲノム構築にあたりオーバーラップしたロングリードを検出するため、グループシードを使用する。 GroupKと呼ばれるこの実装は、小さなオーバーラップまたは両方のリードの高いエラー率によって悪化したオーバーラップを検出するための既存の方法を補完するツールを提供する。この機能により、GroupKは、第3世代のシーケンシングプラットフォームによってシーケンスされたメタゲノムデータにおけるゲノム構築のための賢明な選択を可能にする。これらの集団サンプルは通常不均一な範囲の微生物を含んでいるので、小さなオーバーラップを識別することができることは希少な種のゲノムを再構築するために非常に重要になるだろう。

 

 

インストール

GroupKのこのpythonコードは実験のための古いコードになっている( C ++の実装は今後リリース予定と記載されている)。

  • Python 3 (You may need to modify some script if you want to use Python 2.7)
  • CMake 3.7 or higher
  • GCC 4.8 or higher
  • Biopython library for Python
  • SortedContainers library fot Python
conda install -y sortedcontainers
conda install -y matplotlib
conda install -c bioconda -y sortedcontainers

本体 Github

git clone https://github.com/Strideradu/GroupK.git 
cd GroupK/sa_filter/
g++ -std=c++11 *.cpp -o filter

python GroupK.py -h

# python GroupK.py -h

usage: GroupK.py [-h] [--k1 K1] [--threshold THRESHOLD] [--k2 K2]

                 [--accuracy ACCURACY] [--gap GAP] [--chain CHAIN]

                 [--groupc GROUPC] [--idbase IDBASE] [--ratio RATIO]

                 [--size SIZE] [--large LARGE]

                 input output

 

positional arguments:

  input                 path of input fasta file

  output                path of output file

 

optional arguments:

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

  --k1 K1               kmer size for filtration

  --threshold THRESHOLD

                        count threshold for shared k1

  --k2 K2               kmer size for group

  --accuracy ACCURACY   accuracy of the reads

  --gap GAP             gap rate

  --chain CHAIN         number of kmer in the chain to report

  --groupc GROUPC       the coefficeint c in paper to control chaining

                        threshold, must be non-zero float

  --idbase IDBASE       if the number of identical base meet this threshold

                        the output will always be reported

  --ratio RATIO         ratio of two overlap region threshold

  --size SIZE           group size threshold

  --large LARGE         release threhold for large overlap

usage: GroupK.py [-h] [--k1 K1] [--threshold THRESHOLD] [--k2 K2]

                 [--accuracy ACCURACY] [--gap GAP] [--chain CHAIN]

                 [--groupc GROUPC] [--idbase IDBASE] [--ratio RATIO]

                 [--size SIZE] [--large LARGE]

                 input output

 

positional arguments:

  input                 path of input fasta file

  output                path of output file

 

optional arguments:

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

  --k1 K1               kmer size for filtration

  --threshold THRESHOLD

                        count threshold for shared k1

  --k2 K2               kmer size for group

  --accuracy ACCURACY   accuracy of the reads

  --gap GAP             gap rate

  --chain CHAIN         number of kmer in the chain to report

  --groupc GROUPC       the coefficeint c in paper to control chaining

                        threshold, must be non-zero float

  --idbase IDBASE       if the number of identical base meet this threshold

                        the output will always be reported

  --ratio RATIO         ratio of two overlap region threshold

  --size SIZE           group size threshold

  --large LARGE         release threhold for large overlap

 

 

 

実行方法

FASTA形式のPacbioのシーケンシングリードを指定する。

python GroupK.py --k1 15 --threshold 2 --k2 9\
input.fasta output
  • --k1    kmer size for filtration (default: 15)
  • --threshold    count threshold for shared k1 (default: 2)
  • --k2    kmer size for group (default: 9)
  • --accuracy    accuracy of the reads (default: 0.85)

テスト時はエラーになった。

  

引用

Improving the sensitivity of long read overlap detection using grouped short k-mer matches

Nan Du, Jiao Chen, Yanni Sun
BMC Genomics 2019 20 (Suppl 2) :190

 

 

様々なフォーマットのシーケンスファイルを素早くFASTA形式に変換する any2fasta

 

any2fastaは様々なフォーマットのシーケンスファイルをFASTAフォーマットに変換するperlスクリプトである。他の依存関係はなしにコアのPerlモジュールのみを使用する。非常に高速に実行する。(公開の動機はGithub参照)

 

以下のフォーマットをサポートしている(GIthubより)。

  • Genbank flat file, typically .gb, .gbk, .gbff (starts with LOCUS)
  • EMBL flat file, typically .embl, (starts with ID)
  • GFF with sequence, typically .gff, .gff3 (starts with ##gff)
  • FASTA DNA, typically .fasta, .fa, .fna, .ffn (starts with >)
  • FASTQ DNA, typically .fastq, .fq (starts with @)
  • CLUSTAL alignments, typically .clw, .clu (starts with CLUSTAL or MUSCLE)
  • STOCKHOLM alignments, typically .sth (starts with # STOCKHOLM)
  • GFA assembly graph, typically .gfa (starts with ^[A-Z]\t)

Files may be compressed with:

gzip, typically .gz
bzip2, typically .bz2
zip, typically .zip

 

インストール

macos10.14のminiconda3-4.3.30環境でテストした。

依存

  • any2fasta has no dependencies except Perl 5.10 or higher.

本体 Github

#bioconda
conda install -c bioconda any2fasta

#homebrew
brew install brewsci/bio/any2fasta

#binary
cd /usr/local/bin
wget https://raw.githubusercontent.com/tseemann/any2fasta/master/any2fasta
chmod +x any2fasta

> any2fasta -h

$ any2fasta -h

NAME

  any2fasta 0.4.2

SYNOPSIS

  Convert various sequence formats into FASTA

USAGE

  any2fasta [options] file.{gb,fa,fq,gff,gfa,clw,sth}[.gz,bz2,zip] > output.fasta

OPTIONS

  -h       Print this help

  -v       Print version and exit

  -q       No output while running, only errors

  -n       Replace ambiguous IUPAC letters with 'N'

  -l       Lowercase the sequence

  -u       Uppercase the sequence

HOMEPAGE

  https://github.com/tseemann/any2fasta

END

 

 

テストラン

 入力ファイルと出力ファイルを指定する。

any2fasta input.gbk > output.fasta

 

様々なフォーマットに対応している。NCBI からダウンロードしたgenbankファイルをFASTA形式に変換。

any2fasta GCA_000005845.2_ASM584v2_genomic.gbff.gz > Ecoli.fna 

 

複数ファイルを同時に読み込み、1つのファイルに出力。

any2fasta input1.gbff.gz input2.fna inout3.embl > output.fasta

 

CLUSTALのアラインメントファイルやアセンブリグラフのGFA/FASTGファイルにも対応している。miniasmのアセンブリのGFAファイルをFASTAに変換。

any2fasta assembly.gfa > output.fasta

 

fastqにも対応してますが、処理データ数が多ければseqtk(C lang)やseqkit(Go lang)を使ったほうが早く終わります。

引用

GitHub - tseemann/any2fasta: Convert various sequence formats to FASTA

 

wiki

https://ja.wikipedia.org/wiki/FASTA

 

関連


Torsten Seemannさん(HP)は他にもNGSデータ解析に役立つツールを公開されています(主にスモールゲノム向け)。このブログでもこれまで以下のツールを紹介しました。


 

 

 

アセンブリ配列の16S rRNA相同性からシーケンシングデータの汚染を素早く見積もる ContEst16S

 

 近年、次世代シークエンシング(NGS)と呼ばれる新しいDNAシークエンシング技術の開発により、ゲノムシークエンシングのコストと時間が劇的に減少した。現在、publicデータベースの原核生物ゲノム配列数は約7万に達している(論文執筆時点)。大規模ゲノムデータの使用は、微生物界に関する我々の知識と理解を大いに促進することが示唆されている[ref.1、2]。また、臨床微生物学への応用は、感染症のより良い診断への道を開く[ref.3]。

 NGSの使用が微生物学においてより日常的になるにつれて、汚染を含む配列データの品質保証に関する懸念が高まっている[ref.4–7]。 DNAシーケンシングデータ中の汚染は、生物学的起源(細胞)または試薬もしくは機器中に存在するDNAのいずれかから生じ得る。 NGSは従来のサンガー法よりもはるかに多くの生データ(> 10倍)を生成するため、汚染の可能性が高くなる。汚染は誤った診断につながる可能性があるため、この問題は臨床検査室では特に重要である。そのようなケースを品質管理プロセスとして検出する方法の開発は、日常的な微生物ゲノミクス研究所において最も重要である。

 ドラフトゲノムアセンブリ中の汚染を検出するためのいくつかのアルゴリズムおよびソフトウェアツールが利用可能である。 DeconSeq [ref.8]は、ゲノムまたはメタゲノムアセンブリ中のヒトDNAの検出に特化した、潜在的な汚染物質の事前構築済みデータベースを必要とする。 ProDeGe [ref.9]およびCheckM [ref.10]は、バクテリアおよび古細菌ドメインにわたって高度に保存されているシングルコピータンパク質コード遺伝子を使用する。これらの方法は、publicデータベースのドラフトゲノムアセンブリにおける可能性のある汚染を検出するのに有用である。しかしながら、原則として、それらは汚染を遺伝子移入と区別することができず、それは多くの細菌種においてしばしば起こる[ref.11]。シングルコピータンパク質をコードする遺伝子とは対照的に、rRNA遺伝子は複数のコピーで存在し、水平遺伝子伝播現象を起こしにくいことが知られている[ref.12]。ここで本著者らはContEst16Sと命名した16S rRNA遺伝子配列を用いて原核生物ゲノムアセンブリからの可能性のある生物学的汚染を検出するための新規アルゴリズムを提案する。ここで開発された方法は、publicデータベース内の潜在的に汚染されたゲノムアセンブリを首尾よく同定し、タンパク質をコードする遺伝子に基づく既存のバイオインフォマティクスツールを補完するのに有用であることが証明された。

 


使い方

ContEst16S | Ezbiocloud.net にアクセスする。

f:id:kazumaxneo:20190709004919p:plain

 

Browseボタンから調べたい配列(アセンブルして得たcontig配列)をアップロードする。

f:id:kazumaxneo:20190709004849p:plain

 

バクテリアアーキアを選ぶ。

f:id:kazumaxneo:20190709010556p:plain

Run ContEst16Sボタンを押して実行する。

 

数分で結果が得られる。

コンタミしていることが分かっているデータを使うと、汚染と判定された。

f:id:kazumaxneo:20190709005121p:plain

検出された16S rRNAのAll versus Allペアワイズアラインメント結果。

f:id:kazumaxneo:20190709005130p:plain

 

taxonomic assignment

f:id:kazumaxneo:20190709010804p:plain

 

分子系統樹

f:id:kazumaxneo:20190709005148p:plain

 

引用

ContEst16S: an algorithm that identifies contaminated prokaryotic genomes using 16S RNA gene sequences
Imchang Lee, Mauricio Chalita, Sung-Min Ha, Seong-In Na, Seok-Hwan Yoon, Jongsik Chun

Int J Syst Evol Microbiol. 2017 Jun;67(6):2053-2057.

 

関連