macでインフォマティクス

macでインフォマティクス

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

コンティグの拡張アセンブリを行う PRICE

 

 

低価格のDNAシーケンス技術により、ゲノム、トランスクリプトーム、生態系全体のメタゲノム解析における直接核酸シーケンスの役割は拡大しています。このような大規模なデータセットに対する人間や機械の理解は、配列断片を長く連続した配列ブロック(コンティグ)に合成することで簡略化されますが、アセンブリの分野における進歩のほとんどは、メタゲノムよりもむしろゲノムを単独で対象としています。ここでは、複雑なメタゲノム・データを入力として、特定の核酸種を重点的にアセンブルする戦略であるpaired-read iterative contig extension (PRICE) のためのソフトウェアを紹介する。PRICEは、複雑なメタゲノムから特定の遺伝子、転写産物、ウイルスゲノムを抽出するためのソフトウェアで、カポジ肉腫関連ヘルペスウイルスのBCBL-1株のアセンブリも含まれています。PRICEはオープンソースで、無料でダウンロードできます(derisilab.ucsf.edu/ software/price/ or sourceforge.net/projects/pricedenovo/ )。
 

 

HP (ソースコードが公開されている。マニュアルはリンク先の表をクリックするとダウンロードできる。)

https://derisilab.ucsf.edu/software/price/index.html

PRICE Documentation

https://vcru.wisc.edu/simonlab/bioinformatics/programs/price/PriceDocumentation140408/index.html

 

インストール

ダウンロードしてmakeする。得られるPriceTIに外部依存はないのでそのままパスの通ったディレクトリにコピーする。

ビルド依存

SourceForge

https://sourceforge.net/projects/pricedenovo/files/latest/download

cd PriceSorcexxxx/
make

> ./PriceTI

PRICE Assembler v1.0.1

 

 

These are the command options for the PRICE assembler.

For more details, see the accompanying README.txt file.

Contact: price@derisilab.ucsf.edu

 

Usage: ./PriceTI [args]

 

INPUT FILES:

 accepted formats are fasta (appended .fa, .fasta, .fna, .ffn, .frn),

fastq (.fq, .fastq, or _sequence.txt), or priceq (.pq or .priceq)

  NOTE ABOUT FASTQ ENCODING: multiple encodings are currently used for

fastq quality scores.  The traditional encoding is Phred+33,

                             and PRICE will interpret scores from any

.fq or .fastq file according to that encoding.  The Phred+64

                             encoding has been used extensively by

Illumina, and so it is applied to Illumina's commonly-used

_sequence.txt

                             file append.  Please make sure that your

encoding matches your file append.

INPUT READ FILES:

  NOTE: these flags can be used multiple times in the same command to

include multiple read datasets.

  (default % ID = 90%)

 PAIRED-END FILES (reads are 3p of one another on opposite strands,

i.e. pointing towards one another)

 -fp a b c [d e [f]]: (a,b)input file pair, (c)amplicon insert size

(including read)

                  (d,e,f) are optional; (d)the num. cycles to be

skipped before this file is used;

                  if (f) is provided, then the file will alternate

between being used for (e) cycles and not used for (f) cycles;

                  otherwise, the file will be used for (e) cycles then

will not be used again.

 -fpp a b c d [e f [g]]: (a,b)input file pair, (c)amplicon insert size

(including read), (d)required % identity for match (25-100 allowed)

                  (e,f,g) are optional; (e)the num. cycles to be

skipped before this file is used;

                  if (g) is provided, then the file will alternate

between being used for (f) cycles and not used for (g) cycles;

                  otherwise, the file will be used for (f) cycles then

will not be used again.

 -fs a b [c d [e]]: (a)input paired-end file (alternating entries are

paired-end reads), (b)amplicon insert size (including read)

                  (c,d,e) are optional; (c)the num. cycles to be

skipped before this file is used;

                  if (e) is provided, then the file will alternate

between being used for (d) cycles and not used for (e) cycles;

                  otherwise, the file will be used for (d) cycles then

will not be used again.

 -fsp a b c [d e [f]]: (a)input paired-end file (alternating entries

are paired-end reads), (b)amplicon insert size (including read),

                  (c)required % identity for match (25-100 allowed)

                  (d,e,f) are optional; (d)the num. cycles to be

skipped before this file is used;

                  if (f) is provided, then the file will alternate

between being used for (e) cycles and not used for (f) cycles;

                  otherwise, the file will be used for (e) cycles then

will not be used again.

 MATE-PAIR FILES (reads are 5p of one another on opposite strands,

i.e. pointing away from one another)

 -mp a b c [d e [f]]: like -fp above, but with reads in the opposite

orientation.

 -mpp a b c d [e f [g]]: like -fpp above, but with reads in the

opposite orientation.

 -ms a b [c d [e]]: like -fs above, but with reads in the opposite orientation.

 -msp a b c [d e [f]]: like -fsp above, but with reads in the opposite

orientation.

 FALSE PAIRED-END FILES (unpaired reads are split into paired ends,

with the scores of double-use nuceotides halved)

 -spf a b c [d e [f]]: (a)input file, (b)the length of the 'reads'

that will be taken from each side of the input reads,

                  (c)amplicon insert size (including read)

                  (d,e,f) are optional; (d)the num. cycles to be

skipped before this file is used;

                  if (f) is provided, then the file will alternate

between being used for (e) cycles and not used for (f) cycles;

                  otherwise, the file will be used for (e) cycles then

will not be used again.

 -spfp a b c d [e f [g]]: (a)input file, (b)the length of the 'reads'

that will be taken from each side of the input reads,

                  (c)amplicon insert size (including read),

(d)required % identity for match (25-100 allowed)

                  (e,f,g) are optional; (e)the num. cycles to be

skipped before this file is used;

                  if (g) is provided, then the file will alternate

between being used for (f) cycles and not used for (g) cycles;

                  otherwise, the file will be used for (f) cycles then

will not be used again.

INPUT INITIAL CONTIG FILES:

  NOTE: these flags can be used multiple times in the same command to

include multiple initial contig datasets.

 -icf a b c d: (a)initial contig file, (b)number of addition steps,

(c)number of cycles per step,

               (d)const by which to multiply quality scores

 -picf a b c d e: (a)num of initial contigs from this file, (b)initial

contig file, (c)num addition steps,

                  (d)num cycles per step (e)const by which to multiply

quality scores

 -icfNt/-picfNt: same as -icf/-picf, but if target mode is invoked,

contigs with matches to these input sequences will not necessarily

                 be retained

OUTPUT FILES:

 accepted formats are fasta (.fa or .fasta) or priceq (.pq or .priceq)

 -o a: (a)output file name (.fasta or .priceq)

 -nco a: (a)num. cycles that pass in between output files being

written (default=1)

OTHER PARAMS:

 -nc a: (a)num. of cycles

 -link a: (a)max. number of contigs that are allowed to replace a read

in a mini-assembly (default=2)

 -mol a: (a)minimum overlap length for mini-assembly (default=35)

 -tol a: (a)threshold seq num for scaling overlap for overhang

assemblies (default=20)

 NOTE: -mol and -tol do not affect the parameters for

de-Bruijn-graph-based assembly.

 -mpi a: (a)minimum % identity for mini-assembly (default=85)

 -tpi a: (a)threshold seq num for scaling % ID for mini-assemblies (default=20)

 -MPI, -TPI : same as above, but for meta-assembly (-MPI default=85,

-TPI default=1000)

 NOTE: there is no minimum overlap value for meta-assembly

 -dbmax a: (a) the maximum length sequence that will be fed into de

Bruijn assembly

           (default=100; recommended: max paired-end read length)

 -dbk a: (a) the k-mer size for de Bruijn assembly (default=20; keep

less than the read length)

 -dbms a: (a) the minimum number of sequences to which de Bruijn

assembly will be applied (default=3)

 -r a: (a) alignment score reward for a nucleotide match; should be a

positive integer (default=1)

 -q a: (a) alignment score penalty for a nucleotide mismatch; should

be a negative integer (default=-2)

 -G a: (a) alignment score penalty for opening a gap; should be a

negative integer (default=-5)

 -E a: (a) alignment score penalty for extending a gap; should be a

negative integer (default=-2)

FILTERING READS:

 -rqf a b [c d]: filters pairs of reads if either has an unaccptably

high number of low-quality nucleotides, as defined

                 by the provided quality scores (only applies to files

whose formats include quality score information).

                 (a) the percentage of nucleotides in a read that must

be high-quality; (b) the minimum allowed probability

                 of a nucleotide being correct (must be between 0 and

1, and will usually be a decimal value close to 1);

                 (c) and (d) optionally constrain this filter to use

after (c) cycles have passed, to run for (d) cycles.

                 This flag may be called multiple times to generate

variable behavior across a PRICE run.

 -rnf a [b c]: filters pairs of reads if either has an unaccptably

high number of uncalled nucleotides (Ns or other ambiguous

                 IUPAC codes).  Like -rqf, but will also filter

fasta-format data.  (a) the percentage of nucleotides in a read

                 that must be called; (b) and (c) optionally constrain

this filter to use after (b) cycles have passed, to run

                 for (c) cycles.  This flag may be called multiple

times to generate variable behavior across a PRICE run.

 -maxHp a: filters out a pair of reads if either read has a

homo-polymer track >(a) nucleotides in length.

 -maxDi a: filters out a pair of reads if either read has a repeating

di-nucleotide track >(a) nucleotides in length.

           NOTE: this will also catch mono-nucleotide repeats of the

specified length (a string of A's is also a string

           of AA's), so calling -maxHp in addition to -maxDi is

superfluous unless -maxHp is given a smaller max value.

 -badf a b: prevents reads with a match of at least (b)% identity to a

sequence in file (a) from being mapped to contigs.

 -repmask a b c d e f [g]: uses coverage levels of constructed and/or

input contigs to find repetitive elements and mask them

                           as if they were sequences input using -badf.

                           (a) = cycle number (1-indexed) at which

repeats will be detected.

                           (b) = 's' if repeats will be sought at the

start of the cycle or 'f' if they will be sought at the finish.

                           (c) = the min. number of variance units

above the median that will be counted as high-coverage.

                           (d) = the min. fold increase in coverage

above the median that will be counted as high-coverage.

                           (e) = the min. size in nt for a detected repeat.

                           (f) = reads with a match of at least this %

identity to a repeat will not be mapped to contigs.

                           (g) = an optional output file (.fasta or

.priceq) to which the detected repeats will be written.

 -reset a [b c d...]: re-introduces contigs that were previously not

generating assembly jobs of their own

                      (a) is the one-indexed cycle where the contigs

will be reset.  Same with b, c, d.

                      Any number of args may be added.

FILTERING INITIAL CONTIGS:

 -icbf a b [c]: prevents input sequences with a match of at least (b)%

identity to a sequence in file (a) from being used.

                This filter is optionally not applied to sequences of

length greater than (c) nucleotides.

 -icmHp a [b]: filters out an initial contig if it has a homo-polymer

track >(a) nucleotides in length.

               This filter is optionally not applied to sequences of

length greater than (b) nucleotides.

 -icmDi a [b]: filters out an initial contig if it has a repeating

di-nucleotide track >(a) nucleotides in length.

               NOTE: this will also catch mono-nucleotide repeats of

the specified length (a string of A's is also a string

               of AA's), so calling -icmHp in addition to -icmDi is

superfluous unless -icmHp is given a smaller max value.

               This filter is optionally not applied to sequences of

length greater than (b) nucleotides.

 -icqf a b [c]: filters out an initial contig if it has an unaccptably

high number of low-quality nucleotides, as defined

                by the provided quality scores (only applies to files

whose formats include quality score information).

                (a) the percentage of nucleotides in a read that must

be high-quality; (b) the minimum allowed probability

                of a nucleotide being correct (must be between 0 and

1, and will usually be a decimal value close to 1);

                This filter is optionally not applied to sequences of

length greater than (c) nucleotides.

 -icnf a [b]: filters pairs of reads if either has an unaccptably high

number of uncalled nucleotides (Ns or other ambiguous

              IUPAC codes).  Like -icqf, but will also filter

fasta-format data.  (a) the percentage of nucleotides in a read

              that must be called. This filter is optionally not

applied to sequences of length greater than (c) nucleotides.

FILTERING/PROCESSING ASSEMBLED CONTIGS:

 -lenf a b: filters out contigs shorter than (a) nt at the end of

every cycle, after skipping (b) cycles.

            NOTE: multiple -lenf commands can be entered; for any

cycle, the most recently-initiated filter is used.

            Example: -lenf 50 2 -lenf 300 4 -lenf 200 6 => no filter

for the first two cycles, then a 50nt filter for cycles

                     3 & 4, then a 300nt filter for cycles 5 & 6, then

a 200nt filter for cycles 7 onwards.

 -trim a b [c]: at the end of the (a)th cycle (indexed from 1), trim

off the edges of conigs until reaching the minimum coverage

                level (b), optionally deleting contigs shorter than

(c) after trimming; this flag may be used repeatedly.

 -trimB a b [c]: basal trim; after skipping (a) cycles, trim off the

edges of conigs until reaching the minimum coverage

                 level (b) at the end of EVERY cycle, optionally

deleting contigs shorter than (c) after trimming.

                 -trimB may be called many times, and multiple calls

will interact in the same way as multiple -lenf calls

                 (explained above). A call to -trim will override the

basal trim values for that specified cycle only.

 -trimI a [b]: initial trim; input initial contigs are trimmed before

being used by PRICE to seed assemblies. Contigs are

               trimmed from their outside edges until reaching the

minimum coverage level (a), optionally deleting contigs

               shorter than (b) after trimming. This flag is most

appropriate for .priceq input, can be appropriate for

               .fastq input, and is inappropriate for .fasta input.

It will be equally applied to ALL input contigs.

 -target a b [c d]: limit output contigs to those with matches to

input initial contigs at the end of each cycle.

                    (a) % identity to an input initial contig to count

as a match (ungapped); (b)num cycles to skip

                    before applying this filter.  [c and d are

optional, but must both be provided if either is]

                    After target filtering has begin,

target-filtered/-unfiltered cycles will alternate with (c)

                    filtered cycles followed by (d) unfiltered cycles.

 -targetF a b [c d]: the same as -target, but now matches to all reads

in the input set will be specified, not just

                    the ones that have been introduced up to that

point (this is FullFile mode).

COMPUTATIONAL EFFICIENCY:

 -a x: (x)num threads to use (default=1)

 -mtpf a: (a)max threads per file (default=1)

USER INTERFACE:

 -log a: determines the type of outputmakes the output verbose (lots

of time stamp tags)

         (a) = c: concise stdout (default)

         (a) = n: no stdout

         (a) = v: verbose stdout

 -logf a: (a)the name of an output file for verbose log info to be

written (doesn't change stdout format)

 -, -h, or --help: user interface info.

 

number of cycles MUST be specified with -nc flag

No assembly was run; help message printed instead.

 

 

実行方法

fastaかfastqファイル、もしくはpriceqファイルを指定する。

 

テストラン

パラインフルエンザ4の分離株のゲノムのアセンブルとなっている。

sample job

上のHPからファイルはダウンロードできる。

PriceTI -fpp s_2_1_sequence.txt s_2_2_sequence.txt 300 95 -icf sangerReads.fasta 1 1 5 -nc 30 -dbmax 72 -mol 30 -tol 20 -mpi 80 -target 90 2 1 1 -o practice.fasta

ペアエンドのイルミナFASTQファイルの入力ファイル(Solexaペアエンドリード)を2つ指定している。アンプリコンインサートサイズを300とする。リードをコンティグにアラインする際の配列同一性は95%以上とする。"-icf sangerReads.fasta 1 1 5"の部分はシード配列(sangerリード、454リード、個別のIlluminaリードなど)とパラメータになっている。最初の "1"は、Priceに、入力されたシードの一部を様々なサイクルで追加するのではなく、一度にすべて使用するように指示する。次の "1 "は、別のシードの一部を追加する前に待機するサイクル数を示す。最後の"5"という数字は、品質スコアの定数倍で指定する。"-nc 30"はPRICEに30サイクル実行するよう指示している。"-dbmax 72"はPRICEに、小さな領域のアセンブリにde Bruijn graphサブアセンブラを最大配列長72ヌクレオチドで使うように指示している。"-mol 30"は、配列がマッチする際のオーバーラップの最小値を30ntに設定している。この数値は、-tolフラグで指定された閾値により、カバレッジの増加に応じて自動的にスケールアップされる。多くのアセンブリジョブでは、40ntから50ntの長さで十分である。カバレッジが低いときにこの値を大きくしすぎると、適切なリードがないためにコンティグの成長が止まってしまう可能性がある。"-tol 20"はオーバーラップの最小長をスケーリングするための閾値フラグ。アセンブリジョブのリード数が20を超えると、最小のオーバーラップは上方にスケールされる。"-mpi 80"は、配列が一致するための最小一致率で、80%に設定されている。この値もカバレッジの深さに応じてスケーリングされる。スケーリングを開始する閾値は "-tpi "フラグで指定される。この値を低く設定すると、より緩やかなマッチングが可能になり、値を高く設定すると、マッチングの厳しさが増す。"-target 90 2 1 1"はPRICEを "target "モードで実行することを指定している。これは、最終的なコンティグは入力シード配列を拡張した配列に限定することを意味する。最初の数字 "90 "は、マッチとしてカウントするために必要な、最初のシードのコンティグに対する最小の同一性のパーセンテージである。2番目の数字 "2 "は、PRICEに2サイクル実行させ、アウトボード(ペアエンドの集合によって生成された)コンティグを保持し、出力を最初のシードコンティグとマップまたはコラップしたものに限定するように指示する。次の2つの数字、"1 1"は、その後、アウトボードコンティグの保存とシードコンティグへの制限を交互に行うよう指示している。最後は、"-o practice.fasta"で、PRICEに出力を "practice.fasta "というファイルに格納するよう指示する。

 

実行すると、それぞれのサイクルごとにシードから拡張したfastaファイルが出力される。ファイル名は指定されたprefixとサイクル数を表す"prefix + <INT> + .fasta"となっている。

 

マニュアルより

  • イルミナFASTQを使用する場合、ファイルの末尾は"_sequence.txt "でなければならない。イルミナ以外のFASTQを使用する場合、ファイルの末尾は".fq "または".fastq "でなければならない。
  • FASTAファイルを使用する場合、ファイルの末尾は ".fa" または ".fasta"の必要がある。

 

コメント

渦鞭毛藻類の断片化したプラスチドゲノムを伸ばすために利用されていて気づいたツールです。de novoアセンブリの手法ではうまく組み立てることができない断片化した配列のシード配列からの伸長に利用できそうです。

引用

PRICE: Software for the Targeted Assembly of Components of (Meta) Genomic Sequence Data

J Graham Ruby, Priya Bellare, Joseph L Derisi

G3 (Bethesda). 2013 May 20;3(5):865-80

 

関連

https://kazumaxneo.hatenablog.com/entry/2021/06/17/232020

 

 

 

超高速なfastqの前処理ツール RabbitQCPlus

2023/01/28 追記

 

 シーケンサーデータの品質評価は、ダウンストリームデータ解析において重要な役割を担っている。しかし、既存のツールは、特に圧縮ファイルを扱う場合や、過剰発現解析のような複雑な品質管理操作を行う場合、最適とは言えない効率を達成することが多い。本発表では、最新のマルチコアシステムに対応した超高効率的な品質管理ツールRabbitQCPlusを紹介する。RabbitQCPlusはベクトル化、メモリコピー削減、並列(デ)圧縮、最適化されたデータ構造により、大幅な性能向上を実現している。基本的な品質管理操作を行う場合、最新のアプリケーションと比較して1.1~5.4倍高速になり、必要な計算資源も少なくなる。さらに、RabbitQCPlusはgzip圧縮されたFASTQファイルを処理する場合、他のアプリケーションと比較して少なくとも4倍高速に処理することができる。さらに、280GBのプレーンなFASTQシーケンスデータを処理するのに48コアのサーバーでリードごとの過剰発現解析を有効にすると、他のアプリケーションでは少なくとも22分かかるが、RabbitQCPlusは4分未満で済む。C++ソースコードhttps://github.com/RabbitBio/RabbitQCPlus で入手できる。

 

インストール

ビルド依存

  • gcc 4.8.5 or newer
  • zlib

make時、システムが利用可能なCPUの命令セットは自動認識される。手動で設定することもでき、拡張命令avx256とavx512にも対応している。詳しくはレポジトリを参照。

Github

git clone https://github.com/RabbitBio/RabbitQCPlus.git
cd RabbitQCPlus
make -j4
make install

> RabbitQCPlus -h

RabbitQCPlus

Usage: RabbitQCPlus [OPTIONS]

 

Options:

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

  -i,--inFile1 TEXT           input fastq name 1

  -I,--inFile2 TEXT           input fastq name 2

  -o,--outFile1 TEXT          output fastq name 1

  -O,--outFile2 TEXT          output fastq name 2

  --compressLevel INT         output file compression level (1 - 9), default is 4

  --overWrite                 overwrite out file if already exists.

  --phred64                   input is using phred64 scoring, default is phred33

  --stdin                     input from stdin, or -i /dev/stdin, only for se data or interleaved pe data(which means use --interleavedIn)

  --stdout                    output to stdout, or -o /dev/stdout, only for se data or interleaved pe data(which means use --interleavedOut)

  --notKeepOrder              do not keep order as input when outputing reads (may slightly improve performance if opened), default is off

  -a,--noTrimAdapter          don't trim adapter, default is off

  --decAdaForSe               detect adapter for se data, default is on

  --decAdaForPe               detect adapter for pe data, default is off, tool prefers to use overlap to find adapter

  --printWhatTrimmed          if print what trimmed to *_trimmed_adapters.txt or not, default is not

  --adapterSeq1 TEXT          specify adapter sequence for read1

  --adapterSeq2 TEXT          specify adapter sequence for read2

  --adapterFastaFile TEXT     specify adapter sequences use fasta file

  -c,--correctData            correcting low quality bases using information from overlap, default is off

  -w,--threadNum INT          number of thread used to do QC, including (de)compression for compressed data, default is 8

  -5,--trim5End               do sliding window from 5end to 3end to trim low quality bases, default is off

  -3,--trim3End               do sliding window from 5end to 3end to trim low quality bases, default is off

  --trimFront1 INT            number of bases to trim from front in read1, default is 0

  --trimFront2 INT            number of bases to trim from front in read2, default is 0

  --trimTail1 INT             number of bases to trim from tail in read1, default is 0

  --trimTail2 INT             number of bases to trim from tail in read2, default is 0

  -g,--trimPolyg              do polyg tail trim, default is off

  -x,--trimPolyx              do polyx tail trim, default is off

  -u,--addUmi                 do unique molecular identifier (umi) processing, default is off

  --umiLen INT                umi length if it is in read1/read2, default is 0

  --umiLoc TEXT               specify the location of umi, can be (index1/index2/read1/read2/per_index/per_read), default is 0

  --umiPrefix TEXT            identification to be added in front of umi, default is no prefix

  --umiSkip INT               the number bases to skip if umi exists, default is 0

  --TGS                       process third generation sequencing (TGS) data (only for se data, does not support trimming and will not produce output files), default is off

  -p,--doOverrepresentation   do over-representation sequence analysis, default is off

  -P,--overrepresentationSampling INT

                              do overrepresentation every [] reads, default is 20

  --printORPSeqs              if print overrepresentation sequences to *ORP_sequences.txt or not, default is not

  --noInsertSize              no insert size analysis (only for pe data), default is to do insert size analysis

  --interleavedIn             use interleaved input (only for pe data), default is off

  --interleavedOut            use interleaved output (only for pe data), default is off

  -V,--version                application version

 

 

 

実行方法

ショートリード

#single-end、12スレッド
RabbitQCPlus -w 12 -i in1.fastq.gz -o p1.fastq.gz

#paired-end、20スレッド
RabbitQCPlus -w 40 -i in1.fastq.gz -I in2.fastq.gz -o p1.fastq.gz -O p2.fastq.gz

 

 

ロングリード

RabbitQCPlus -w 4 -i in.fastq --TGS

 

オプションを指定しない場合もhtmlレポートも出力される。レポジトリの例。

 

コメント

ファイルサイズ12GBの生のペアエンドのfastqファイルの前処理にかかった時間は、

64 CPU - 12.5秒

40 CPU - 15.0秒

20 CPU - 31.7秒

という結果でした(TR3990x、256GBメモリ、PCIe4接続SSD環境、n=1)。

このストレージ上で同じファイルをコピーして書き出すとおよそ10.5秒かかりました(ubuntuのcpコマンド使用)。12.5秒というタイムは理論限界値に近い値が出ていると言えるんじゃないかと思います(注;本SSDのシーケンシャル読み書きのスペック値は6GB/s)。avx256やavx512拡張命令に対応した最近のメニーコアCPUとPCIe5以上の複数レーンに対応した高速なストレージを組み合わせると、より時間短縮できるかもしれません。

引用

RabbitQCPlus: More Efficient Quality Control for Sequencing Data

2022 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)

DOI Bookmark: 10.1109/BIBM55620.2022.9995332

 

関連


 

プロトコルを共有する Bio-protocol Exchange

 

Bio-protocol Exchangeは、Bio-protocolジャーナルの姉妹プラットフォームサイト。研究者がプロトコルを議論したり、実験手順を共有する場として公開されている。実験プロトコルだけのプレプリントサーバーとしても機能している。

 

 

editorial team

https://bio-protocol.org/exchange/editors

オーサーガイドライン

Bio-protocol Exchange, a platform for researchers to find, share and discuss life science protocols

 

HP

Bio-protocol Exchange, a platform for researchers to find, share and discuss life science protocols

 

生物学実験のプロトコルについて質問・共有、意見交換したり、バイオインフォマティクスの実験プロトコルを検索することができる。

 

上のメニューからはBio-protocolジャーナルにもリンクしている。

 

ここではprotocolのページからmakerと検索してみる。補足すると、このMAKERは真核生物のアノテーションパイプラインで様々なツールを統合している。 ab initio遺伝子予測、エビデンスベースの遺伝子予測、EVMによる結果の統合、ゲノムブラウザ使用など柔軟に対応する。しかしその分ベストプラクティスなワークフローを見つけるのは簡単ではなく、スモールラボなどで情報交換が難しい人にとって、適切な使い方を習得するまでの労力が大きくなってしまっている。

 

4816ヒットした。Peer-reviwedはBio-protocolジャーナルから出版された論文へのヒット。右のConsiceがプロトコルや質問の検索結果。



検索結果は右上のsort byからソートできる。

 

 

1つ開いてみる。上に論文へのリンクが張ってあり、すでに論文になっている研究のプロトコルであることが分かる。DOIもついている(引用可能であり、削除することはできない)。

質問もできるようになっている。”明確で、具体的、かつ簡潔な質問をお願いします”と注意書きもある。

 

あまり情報が見つからない(バイオ)インフォマティクスの実験プロトコルを探すためにも利用できると思います。アクセスしてみてください。

 

 

3'UTRのアノテーションを行う peaks2utr

 

 非モデル生物のアノテーションは未解決の問題であり、特に非翻訳領域(UTR)の検出が重要である。UTRの正確なアノテーションはトランスクリプトーム解析において各遺伝子の発現を正確に把握するために非常に重要であるが、アノテーションパイプラインではほとんど見落とされているのが現状である。本ツールは、10x ChromiumのようなシングルセルテクノロジーのUTRエンリッチメントを利用して、与えられたカノニカルアノテーションに対して3' UTRを正確にアノテーションする、使いやすいPythonコマンドラインツールpeaks2utrを紹介する。Python 3 (>= 3.8) で実装されており、PyPI (https://pypi.org/project/peaks2utr) および GitHub (https://github.com/haessar/peaks2utr) で入手できる。

 

GIthubより

peaks2utrは、BAM形式のアライメント済みリードと、GFFまたはGTF形式の標準的なアノテーションのセットに対して、3'非翻訳領域(UTR)をアノテーションする。peaks2utrでは、MACS2 (https://pypi.org/project/MACS2/) を使ってBAMファイル内の有意なリードカバレッジの幅広い「ピーク」を呼び出し、一連の基準に合格したピークを用いて新規3' UTRのアノテーションを付ける。これは、3'または5'UTRの遠位端にシグナルが集中する10x ChromiumランのようなBAMファイルを優先する。ソフトクリップ塩基や所定の長さのpolyA-tailを含むリードを検出し、それらの末端塩基を「切断点」として集計する。peaks2utrは、入力GFFファイル中の既存の3'UTRアノテーションを拡張、上書き、無視するように調整することができる。

インストール

Github

#pip
pip install peaks2utr

#install from source
git clone https://github.com/haessar/peaks2utr.git
cd peaks2utr

> peaks2utr -h

usage: peaks2utr [-h] [--max-distance MAX_DISTANCE] [--override-utr]

[--extend-utr] [--five-prime-ext FIVE_PRIME_EXT] [--skip-soft-clip]

[--min-pileups MIN_PILEUPS] [--min-poly-tail MIN_POLY_TAIL] [-p

PROCESSORS] [-f] [-o OUTPUT] [--keep-cache] GFF_IN BAM_IN

 

Use MACS2 to build forward and reverse peaks files for given .bam

file. Iterate peaks through set of criteria to determine UTR

viability, before annotating in .gff file.

 

positional arguments:

  GFF_IN                input 'canonical' annotations file in gff or gtf format.

  BAM_IN                input reads file in bam format.

 

optional arguments:

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

  --max-distance MAX_DISTANCE

                        maximum distance in bases that UTR can be from

a transcript.

  --override-utr        ignore already annotated 3' UTRs in criteria.

  --extend-utr          extend previously existing 3' UTR annotations

where possible.

  --five-prime-ext FIVE_PRIME_EXT

                        a peak within this many bases of a gene's

5'-end should be assumed to belong to it.

  --skip-soft-clip      skip the resource-intensive logic to pileup

soft-clipped read edges.

  --min-pileups MIN_PILEUPS

                        Minimum number of piled-up mapped reads for UTR cut-off.

  --min-poly-tail MIN_POLY_TAIL

                        Minimum length of poly-A/T tail considered in

soft-clipped reads.

  -p PROCESSORS, --processors PROCESSORS

                        How many processor cores to use.

  -f, -force, --force   Overwrite outputs if they exist.

  -o OUTPUT, --output OUTPUT

                        output filename.

  --keep-cache          Keep cached files on run completion.

 

 

 

テストラン

wget https://github.com/haessar/peaks2utr/raw/master/demo/Tb927_01_v5.1.gff
wget https://github.com/haessar/peaks2utr/raw/master/demo/Tb927_01_v5.1.slice.bam
peaks2utr Tb927_01_v5.1.gff Tb927_01_v5.1.slice.bam

出力例

 

8CPU指定

peaks2utr -p 8 -o output.gff3 canonical.gff reads.bam

 

注;RNA seqはUTRの末端が読まれるような特別なプロトコルで実行されている必要がある。そうでないRNA seqからのbamを使うと不完全なUTR予測結果になる。

引用

peaks2utr: a robust Python tool for the annotation of 3' UTRs

William Haese-Hill,  Kathryn Crouch,  Thomas D Otto

bioRxiv, Posted May 26, 2022

 

参考

Ensembl insights: How are UTRs annotated?

https://www.ensembl.info/2018/08/17/ensembl-insights-how-are-utrs-annotated/

 

アノテーションパイプライン MAKER

2008年の論文

 移植可能で容易に設定可能なゲノムアノテーションパイプラインであるMAKERを開発した。MAKERの目的は、研究者が独立して真核生物ゲノムのアノテーションを行い、ゲノムデータベースを作成することである。MAKERはリピートを識別し、ESTやタンパク質をゲノムにアラインさせ、ab initio遺伝子予測を行い、これらのデータを証拠に基づく品質指標を持つ遺伝子アノテーションに自動的に生成する。また、MAKERは簡単にトレーニングすることができる。予備実験の結果をもとに遺伝子予測アルゴリズムを自動的に再トレーニングし、次回以降、より質の高い遺伝子モデルを生成することができる。MAKERの入力は最小限であり、その出力はGMODデータベースの作成に使用することができる。MAKERの出力はApollo Genomeブラウザで見ることができる。この機能により、データベースのオーバーヘッドなしに、個々のコンティグやBACのアノテーション、表示、編集を簡単に行うことができるようになる。MAKERは、プラナリアのSchmidtea mediterraneaのゲノムをアノテーションし、新しいゲノムデータベースSmedGDを作成するために使用された。また、MAKERの性能を他のアノテーションパイプラインと比較した。その結果、MAKERはゲノム配列をコミュニティがアクセス可能なゲノムデータベースに変換するための簡単で効果的な手段を提供することが実証された。MAKERは、広範なバイオインフォマティクスリソースが容易に利用できない新興のモデル生物ゲノムプロジェクトに特に有用である。

チュートリアルより

MAKERは、バイオインフォマティクスの経験が少ない小規模な研究グループでも使用できるように設計された、使いやすいゲノムアノテーションパイプラインである。MAKERはスケーラブルに設計されており、大規模なシーケンシングセンターを含むプロジェクトにも適している。

MAKER はアノテーションパイプラインであり、遺伝子予測ツールではない。MAKERは遺伝子を予測するのではなく、既存のソフトウェアツール(一部は遺伝子予測ツール)を活用し、エビデンスアラインメントに基づいて最適な遺伝子モデルを生成する。(注;MAKER自身は遺伝子予測ツールではないが、内部でab initio遺伝子予測ツールを訓練し、遺伝子予測すること自体は可能。tutorialの後半参照)

 

What does MAKER do?

  • Identifies and masks out repeat elements
  • Aligns ESTs to the genome
  • Aligns proteins to the genome
  • Produces ab initio gene predictions
  • Synthesizes these data into final annotations
  • Produces evidence-based quality values for downstream annotation management

 

HP

https://www.yandell-lab.org/software/maker.html

 

MAKER Tutorial for WGS Assembly and Annotation Winter School 2018 

http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/MAKER_Tutorial_for_WGS_Assembly_and_Annotation_Winter_School_2018

 

インストール

Github

#conda(link) 
mamba install -c "bioconda/label/cf201901" -c conda-forge -c bioconda maker
=>condaでインストールした場合もRepbaseからrepeatmakerのRepbaseのrepeatmakerライブラリをダウンロードして配置する必要がある*1

#Docker image
#ここでは公開されているdocker imageを利用する(link)1-5の手順を踏む
#1 base image(link)のpull。このイメージ自体は4つの親イメージから順に作成されている。ビルドしても良い。
docker pull chrishah/premaker-plus:18-0d9787e

#2 作業ディレクトリでmakerのレポジトリをpullする
cd <change>/<to>/<working_dir>
git clone https://github.com/Yandell-Lab/maker

#3 RepbaseからrepeatmakerのRepbaseのrepeatmakerライブラリをとってくる(link)。最新の全データダウンロードには登録だけではだめで、所属機関のライセンスが必要(有料。昔から?)。ないのでここでは古めだがこれを使用(とても助かります)
=> tar ballを解凍し、後でRepeatMasker/Libraries/に配置する

#4 Dockerfile作成とビルド(#1のbase imageから)
echo -e "FROM chrishah/premaker-plus:18-0d9787e" > Dockerfile-maker-plus
docker build --network=host -t maker-plus:2.31 --file Dockerfile-maker-plus .

#5 出来たimageをランしてMaker(#2でcloneしたレポジトリ)をビルドする
docker run -itv $PWD:/home -w /home maker-plus:2.31
$ cd /home/maker/src
$ perl Build.PL
$ ./Build install
$ cd ../bin
$ ./maker #必要ならmaker/bin/にパスを通す

> maker -h

MAKER version 3.01.03

 

Usage:

 

     maker [options] <maker_opts> <maker_bopts> <maker_exe>

 

 

Description:

 

     MAKER is a program that produces gene annotations in GFF3 format using

     evidence such as EST alignments and protein homology. MAKER can be used to

     produce gene annotations for new genomes as well as update annotations

     from existing genome databases.

 

     The three input arguments are control files that specify how MAKER should

     behave. All options for MAKER should be set in the control files, but a

     few can also be set on the command line. Command line options provide a

     convenient machanism to override commonly altered control file values.

     MAKER will automatically search for the control files in the current

     working directory if they are not specified on the command line.

 

     Input files listed in the control options files must be in fasta format

     unless otherwise specified. Please see MAKER documentation to learn more

     about control file  configuration.  MAKER will automatically try and

     locate the user control files in the current working directory if these

     arguments are not supplied when initializing MAKER.

 

     It is important to note that MAKER does not try and recalculated data that

     it has already calculated.  For example, if you run an analysis twice on

     the same dataset you will notice that MAKER does not rerun any of the

     BLAST analyses, but instead uses the blast analyses stored from the

     previous run. To force MAKER to rerun all analyses, use the -f flag.

 

     MAKER also supports parallelization via MPI on computer clusters. Just

     launch MAKER via mpiexec (i.e. mpiexec -n 40 maker). MPI support must be

     configured during the MAKER installation process for this to work though

     

 

Options:

 

     -genome|g <file>    Overrides the genome file path in the control files

 

     -RM_off|R           Turns all repeat masking options off.

 

     -datastore/         Forcably turn on/off MAKER's two deep directory

      nodatastore        structure for output.  Always on by default.

 

     -old_struct         Use the old directory styles (MAKER 2.26 and lower)

 

     -base    <string>   Set the base name MAKER uses to save output files.

                         MAKER uses the input genome file name by default.

 

     -tries|t <integer>  Run contigs up to the specified number of tries.

 

     -cpus|c  <integer>  Tells how many cpus to use for BLAST analysis.

                         Note: this is for BLAST and not for MPI!

 

     -force|f            Forces MAKER to delete old files before running again.

             This will require all blast analyses to be rerun.

 

     -again|a            recaculate all annotations and output files even if no

             settings have changed. Does not delete old analyses.

 

     -quiet|q            Regular quiet. Only a handlful of status messages.

 

     -qq                 Even more quiet. There are no status messages.

 

     -dsindex            Quickly generate datastore index file. Note that this

                         will not check if run settings have changed on contigs

 

     -nolock             Turn off file locks. May be usful on some file systems,

                         but can cause race conditions if running in parallel.

 

     -TMP                Specify temporary directory to use.

 

     -CTL                Generate empty control files in the current directory.

 

     -OPTS               Generates just the maker_opts.ctl file.

 

     -BOPTS              Generates just the maker_bopts.ctl file.

 

     -EXE                Generates just the maker_exe.ctl file.

 

     -MWAS    <option>   Easy way to control mwas_server for web-based GUI

 

                              options:  STOP

                                        START

                                        RESTART

 

     -version            Prints the MAKER version.

 

     -help|?             Prints this usage statement.

 

fasta_merge 

Synopsis:

 

fasta_merge -d maker_datastore_index.log

fasta_merge -o genome.all -i <fasta1> <fasta2> ...

 

Descriptions:

 

This script will take a MAKER datastore index log file, extract all

the relevant fasta files and create fasta files with relevant

categories of sequence (i.e. transcript, protein, GeneMark protien,

etc.).  For this to work properly you need to be in the same directory

as the datastore index.

 

Options:

 

  -d The location of the MAKER datastore index log.

  -o Alternate base name for the output files.

  -i A optional list of files to process along with or instead of the

     datastore.

 

> gff3_merge 

Synopsis:

 

gff3_merge -d maker_datastore_index.log

gff3_merge -o genome.all.gff <gff3_file1> <gff3_file2> ...

 

Descriptions:

 

This script will take a MAKER datastore index log file, extract all

the relevant GFF3 files and combined GFF3 file.  The script can also

combine other correctly formated GFF3 files.  For this to work

properly you need to be in the same directory as the datastore index.

 

Options:

 

  -d The location of the MAKER datastore index log file.

  -o Alternate base name for the output files.

  -s Use STDOUT for output.

  -g Only write MAKER gene models to the file, and ignore evidence.

  -n Do not print fasta sequence in footer

  -l Merge legacy annotation sets (ignores already having seen

     features more than once for the same contig)

 

 

makerライセンスに関する注意事項

 

 

テストラン

レポジトリのdata/にテスト用のデータセットが用意されている。

cd maker/data/

maker/data/

1,まずconfigファイルのtempleteを作成する。MAKER は現在のパスにあるconfigファイルを(明示的に指示することなく)自動で探して実行するようになっている。そのため、各ゲノム固有のディレクトリで MAKER を実行することが推奨されている。

maker -CTL

4つのファイルが出来る。

マニュアルの説明

  • maker_exe.ctl - 基礎となる実行ファイルのパス情報が含まれています。
  • maker_bopt.ctl - BLASTとExonerateのフィルタリングの統計情報が含まれています。
  • maker_opt.ctl - 入力ゲノムファイルの場所を含む、MAKERに関するその他のすべての情報が含まれています。

このあと、マニュアルではnanoエディタ(apt install nano)を使ってconfigファイルを編集している。

 

maker_exe.ctl(実行ファイルのパス情報を含む)を開いてみる。

上のような手順で正しくインストール手順に従った場合、すべての実行ファイルのパスが表示されているはず。編集する場合は"="の両側にはスペースを入れてはいけない。

 

maker_opt.ctl(入力ファイルのパスを含むその他のすべての情報)を開く。

maker_opt.ctlの6つの行を以下のように修正する。

genome=dpp_contig.fasta
est=dpp_est.fasta
protein=dpp_protein.fasta
est2genome=1

protein2genome=1

 

これで準備ができた。

 

3,編集が終わったらMAKERをランする。(*1)

maker -q -base annotation
  •  -q     Regular quiet. Only a handlful of status messages.
  • -qq    Even more quiet. There are no status messages.
  • -base     Set the base name MAKER uses to save output files. MAKER uses the input genome file name by default.
  •  -R    Turns all repeat masking options off.
  • -f       Forces MAKER to delete old files before running again.  This will require all blast analyses to be rerun.

  • -a      recaculate all annotations and output files even if no settings have changed. Does not delete old analyses.

     

...

 

出力

dpp_contig.maker.output/ #入力ゲノムファイル名に基づいている

 

maker_opts.log, maker_exe.log, maker_bopts.log ファイルは、今回のMAKERの実行に使用された制御ファイルのログ。mpi_blastdbディレクトリには、入力されたEST、タンパク質、リピートデータベースから作成されたFASTAインデックスとBLASTデータベースファイルが入っている。

dpp_contig_datastoreディレクトリには、ゲノムFastaファイルから個々のコンティグに対する最終的なMAKER出力を格納したサブフォルダ群があります。

 

dpp_contig_master_datastore_index.logがログファイルになっている。

> cat dpp_contig.maker.output/dpp_contig_master_datastore_index.log

サンプルファイルにはコンティグが1つしかなかったため、1つのコンティグを記述したエントリしかない。写真は、コンティグcontig-dpp-500-500がSTARTし、その後、何事もなくFINISHEDしたことを示している。さらにこのコンティグの結果がdpp_contig_datastore/05/1F/contig-dpp-500-500/に保存されたことも示している。(マニュアルより)。ほかにも以下の種類がある。

FAILED - このコンティグの実行に失敗
RETRY - MAKERが失敗したコンティグを再試行した。
SKIPPED_SMALL - コンティグが短すぎてアノテーションができない(最小長は maker_opt.ctl で指定できる)。

DIED_SKIPPED_PERMANENT - MAKERが再試行でも失敗したコンティグ(コンティグの再試行回数は maker_opt.ctl で指定できる)。

数千から数十万のコンティグが含まれる場合、ネットワーク経由のアクセスでもパフォーマンスが低下する可能性がある。MAKERは「ベース」から始まるネストしたサブディレクトリの階層を作成し、与えられたコンティグの結果を、ネストした数千のディレクトリのデータストアに配置し、これを回避する。master_datastore_index.logはンティグの出力がどこに格納されているかを特定するためにも不可欠なものになる(マニュアルより))。

 

コンティグ1つの結果を見てみる。

dpp_contig.maker.output/dpp_contig_datastore/05/1F/contig-dpp-500-500/

contig-dpp-500-500.gffは、GFF3形式のアノテーションファイル。contig-dpp-500-500.maker.transcripts.fasta と contig-dpp-500-500.maker.proteins.fasta ファイルには、遺伝子アノテーションの転写産物とタンパク質の配列が保存されている。theVoid.contig-dpp-500-500ディレクトリには、MAKERが実行するすべてのプログラム(Blast、SNAP、RepeatMaskerなど)からの個別の出力ファイルが保存されている。

 

 

4、各コンティグの結果のマージ

cd dpp_contig.maker.output/
fasta_merge -d dpp_contig_master_datastore_index.log
=> dpp_contig.all.maker.transcripts.fastaができる

gff3_merge -d dpp_contig_master_datastore_index.log
=> dpp_contig.all.gffができる

(注;テストデータにはコンティグが1つしかない)

 

こちらの論文ではMAKERを使って2段階のアノテーションを行うための設定が書かれています。参考になります。

An improved genome assembly uncovers prolific tandem repeats in Atlantic cod | BMC Genomics | Full Text

 

注;MAKERはMPIに対応していて、並列計算することで速度を大幅にアップすることできます。あるサイズの領域に切り分けて、それぞれの領域ごとに独立して計算できるからだと思われます。言い換えれば、単一ノードで実行するとかなり膨大な時間がかかるということです。ある程度大きなゲノムを扱う場合(例えば100Mb以上)、注意してください。

引用

MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes
Brandi L Cantarel, Ian Korf, Sofia M C Robb, Genis Parra, Eric Ross, Barry Moore, Carson Holt, Alejandro Sánchez Alvarado, Mark Yandell

Genome Res. 2008 Jan;18(1):188-96

MAKER-P: a tool kit for the rapid creation, management, and quality control of plant genome annotations
Michael S Campbell, MeiYee Law, Carson Holt, Joshua C Stein, Gaurav D Moghe, David E Hufnagel, Jikai Lei, Rujira Achawanantakun, Dian Jiao, Carolyn J Lawrence, Doreen Ware, Shin-Han Shiu, Kevin L Childs, Yanni Sun, Ning Jiang, Mark Yandell

Plant Physiol. 2014 Feb;164(2):513-24. doi: 10.1104/pp.113.230144. Epub 2013 Dec 4.

 

Genome Annotation and Curation Using MAKER and MAKER-P
Michael S. Campbell, Carson Holt, Barry Moore, Mark Yandell
Curr Protoc Bioinformatics. 2014 Dec 12;48:4.11.1-4.11.39

PMCアーカイブlink

 

参考

*1

Repbaseのライブラリがないと怒られる。それでも強制的にランするにはリピートマスク処理をOFFにする(-R)か、maker_opts.ctlのmodel_orgを写真のように空白にする。

bioconda-recipes:issue

ERROR: Could not determine if RepBase is installed · Issue #16501 · bioconda/bioconda-recipes · GitHub

 

もしくはrmlib=にfasta形式のリピートライブラリファイルを指定する(現在のバージョンのRepeatMakerはfasta形式のリピートファイルを認識する)。例えばRepatModelerを使えば種固有の性質があるリピートをDe novo探索し、fasta形式で出力する、リピートライブラリがない新しいゲノムではこれを使用できる(紹介)。その場合、"model_org="を空白にし、"rmlib="にRepatModelerで予測したリピートのコンセンサス配列のfastaファイルを指定する。

 

関連


 

 

 

 

cDNA配列をゲノムにアラインメントする GMAP

 

2016年の論文より

RNA-SeqやDNA-SeqのデータセットをゲノムにアライメントするためのプログラムGMAPとGSNAPは、生物学的手法の進歩に伴い、より長いリード、より大量のデータ、新しいタイプの生物学的アッセイを扱うために進化してきた。ゲノム表現では、SIMD(single-instruction multiple-data)命令を用いて配列を比較できる線形ゲノム、SIMD命令を用いて高速アクセスできる圧縮ゲノムハッシュテーブル、40億bp以上の大規模ゲノムへの対応、高速アクセスのための新規データ構造を持つ拡張サフィックスアレイ(ESA)などに改良されてきた。アルゴリズムの改良では、接尾辞配列を用いた greedy match-and-extend アルゴリズム、ゲノムハッシュテーブルを用いたセグメントチェイニング、セグメントハッシュテーブルを用いたdiagonalization、SIMD命令を用いた核酸レベルの動的計画法(Fループ計算を不要にする)などが行われた。プログラムの機能強化には、indelの位置の標準化、あいまいなスプライシングの処理、重複するペアエンドリードのクリッピングとマージ、環状染色体やスキャフォールドへのアラインメントが含まれる。これらのプログラムは、gmapR や HTSeqGenie などの R/Bioconductor パッケージに統合され、パイプラインで使用できるようになり、これらのパイプラインは多くの生命現象の発見を促進している。

 

インストール

Github

git clone https://github.com/juliangehring/GMAP-GSNAP.git
cd GMAP-GSNAP/
./configure
make
make check   (optional)
make install

#conda
mamba install -c bioconda gmap

> gmap

$ gmap

Note: gmap.avx2 does not exist.  For faster speed, may want to compile package on an AVX2 machine GMAP version 2021-08-25 called with args: gmap.sse42 Need to specify the -d, -g, -1, -2, or --cmdline flag Usage: gmap [OPTIONS...] <FASTA files...>, or

       cat <FASTA files...> | gmap [OPTIONS...]

 

Input options (must include -d or -g)

  -D, --dir=directory            Genome directory.  Default (as specified by --with-gmapdb to the configure program) is mambaforge/envs/gemoma/share

  -d, --db=STRING                Genome database.  If argument is '?' (with the quotes), this command lists available databases.

 

  -k, --kmer=INT                 kmer size to use in genome database (allowed values: 16 or less).

                                   If not specified, the program will find the highest available

                                   kmer size in the genome database

  --sampling=INT                 Sampling to use in genome database. If not specified, the program

                                   will find the smallest available sampling value in the genome database within selected k-mer size

  -g, --gseg=filename            User-supplied genomic segment

  -1, --selfalign                Align one sequence against itself in FASTA format via stdin

                                   (Useful for getting protein translation of a nucleotide sequence)

  -2, --pairalign                Align two sequences in FASTA format via stdin, first one being

                                   genomic and second one being cDNA

  --cmdline=STRING,STRING        Align these two sequences provided on the command line,

                                   first one being genomic and second one being cDNA

  -q, --part=INT/INT             Process only the i-th out of every n sequences

                                   e.g., 0/100 or 99/100 (useful for distributing jobs

                                   to a computer farm).

  --input-buffer-size=INT        Size of input buffer (program reads this many sequences

                                   at a time for efficiency) (default 1000)

 

Computation options

  -B, --batch=INT                Batch mode (default = 2)

                                 Mode     Positions       Genome

                                   0      mmap            mmap

                                   1      mmap & preload  mmap

                      (default)    2      mmap & preload  mmap & preload

                                   3      allocate        mmap & preload

                                   4      allocate        allocate

                                   5      allocate        allocate (same as 4)

                           Note: For a single sequence, all data structures use mmap If mmap not available and allocate not chosen, then will use fileio (very slow)

  --use-shared-memory=INT        If 1, then allocated memory is shared among all processes on this node

                                   If 0 (default), then each process has private allocated memory

  --nosplicing                   Turns off splicing (useful for aligning genomic sequences

                                   onto a genome)

  --max-deletionlength=INT       Max length for a deletion (default 100).  Above this size,

                                   a genomic gap will be considered an intron rather than a deletion.

                                   If the genomic gap is less than --max-deletionlength and greater

                                   than --min-intronlength, a known splice site or splice site probabilities

                                   of 0.80 on both sides will be reported as an intron.

  --min-intronlength=INT         Min length for one internal intron (default 9).  Below this size,

                                   a genomic gap will be considered a deletion rather than an intron.

                                   If the genomic gap is less than

--max-deletionlength and greater

                                   than --min-intronlength, a known splice site or splice site probabilities

                                   of 0.80 on both sides will be

reported as an intron.

  --max-intronlength-middle=INT  Max length for one internal intron (default 500000).  Note: for backward compatibility, the -K or

--intronlength flag will set both --max-intronlength-middle and

--max-intronlength-ends. Also see --split-large-introns below.

  --max-intronlength-ends=INT    Max length for first or last intron (default 10000).  Note: for backward  compatibility, the -K or

--intronlength flag will set both  --max-intronlength-middle and --max-intronlength-ends.

  --split-large-introns          Sometimes GMAP will exceed the value for --max-intronlength-middle,  if it finds a good single alignment.  However, you can force GMAP to split such alignments by using this flag

  --trim-end-exons=INT           Trim end exons with fewer than given number of matches

 (in nt, default 12)

  -w, --localsplicedist=INT      Max length for known splice sites at ends of sequence

                                   (default 2000000)

  -L, --totallength=INT          Max total intron length (default 2400000)

  -x, --chimera-margin=INT       Amount of unaligned sequence that triggers

                                   search for the remaining sequence (default 30). Enables alignment of chimeric reads, and may help with some non-chimeric reads.  To turn off, set to zero.

  --no-chimeras                  Turns off finding of chimeras.  Same effect as --chimera-margin=0

  -t, --nthreads=INT             Number of worker threads

  -c, --chrsubset=string         Limit search to given chromosome

  --strand=STRING                Genome strand to try aligning to (plus, minus, or both default)

  -z, --direction=STRING         cDNA direction (sense_force, antisense_force,

                                   sense_filter, antisense_filter,or auto (default))

  --canonical-mode=INT           Reward for canonical and semi-canonical introns

                                   0=low reward, 1=high reward (default), 2=low reward for

                                   high-identity sequences and high reward otherwise

  --cross-species                Use a more sensitive search for canonical splicing, which helps especially

                                   for cross-species alignments and other difficult cases

  --allow-close-indels=INT       Allow an insertion and deletion close to each other

                                   (0=no, 1=yes (default), 2=only for

high-quality alignments)

  --microexon-spliceprob=FLOAT   Allow microexons only if one of the splice site probabilities is greater than this value (default 0.95)

  --indel-open                   In dynamic programming, opening penalty for indel

  --indel-extend                 In dynamic programming, extension penalty for indel

                                   Values for --indel-open and

--indel-extend should be in [-127,-1].

                                   If value is < -127, then will use -127 instead.

                                   If --indel-open and --indel-extend are not specified, values are chose adaptively, based on the differences between the query and reference

  --cmetdir=STRING               Directory for methylcytosine index files (created using cmetindex)

                                   (default is location of genome index files specified using -D, -V, and -d)

  --atoidir=STRING               Directory for A-to-I RNA editing index files (created using atoiindex)

                                   (default is location of genome index files specified using -D, -V, and -d)

  --mode=STRING                  Alignment mode: standard (default), cmet-stranded, cmet-nonstranded,

                                    atoi-stranded, atoi-nonstranded, ttoc-stranded, or ttoc-nonstranded.

                                    Non-standard modes requires you to have previously run the cmetindex

                                    or atoiindex programs (which also cover the ttoc modes) on the genome

  -p, --prunelevel               Pruning level: 0=no pruning (default), 1=poor seqs,

                                   2=repetitive seqs, 3=poor and repetitive

 

Output types

  -S, --summary                  Show summary of alignments only

  -A, --align                    Show alignments

  -3, --continuous               Show alignment in three continuous lines

  -4, --continuous-by-exon       Show alignment in three lines per exon

  -E, --exons=STRING             Print exons ("cdna" or "genomic")

                                   Will also print introns with

"cdna+introns" or

                                   "genomic+introns"

  -P, --protein_dna              Print protein sequence (cDNA)

  -Q, --protein_gen              Print protein sequence (genomic)

  -f, --format=INT               Other format for output (also note the -A and -S options

                                   and other options listed under Output types):

                                   mask_introns,

                                   mask_utr_introns,

                                   psl (or 1) = PSL (BLAT) format,

                                   gff3_gene (or 2) = GFF3 gene format,

                                   gff3_match_cdna (or 3) = GFF3 cDNA_match format,

                                   gff3_match_est (or 4) = GFF3 EST_match format,

                                   splicesites (or 6) = splicesites output (for GSNAP splicing file),

                                   introns = introns output (for GSNAP splicing file),

                                   map_exons (or 7) = IIT FASTA exon map format,

                                   map_ranges (or 8) = IIT FASTA range map format,

                                   coords (or 9) = coords in table format,

                                   sampe = SAM format (setting paired_read bit in flag),

                                   samse = SAM format (without setting paired_read bit),

                                   bedpe = indels and gaps in BEDPE format Output options

  -n, --npaths=INT               Maximum number of paths to show (default 5).  If set to 1, GMAP

                                   will not report chimeric alignments, since those imply

                                   two paths.  If you want a single alignment plus chimeric

                                   alignments, then set this to be 0.

  --suboptimal-score=FLOAT       Report only paths whose score is within this value of the

                                   best path.

                                 If specified between 0.0 and 1.0, then treated as a fraction

                                   of the score of the best alignment (matches minus penalties for

                                   mismatches and indels).  Otherwise, treated as an integer

                                   number to be subtracted from the score of the best alignment.

                                   Default value is 0.50.

  -O, --ordered                  Print output in same order as input (relevant

                                   only if there is more than one worker thread)

  -5, --md5                      Print MD5 checksum for each query sequence

  -o, --chimera-overlap          Overlap to show, if any, at chimera breakpoint

  --failsonly                    Print only failed alignments, those with no results

  --nofails                      Exclude printing of failed alignments

 

  -V, --snpsdir=STRING           Directory for SNPs index files (created using snpindex) (default is location of genome index files specified using -D and -d)

   -v, --use-snps=STRING          Use database containing known SNPs (in <STRING>.iit, built

                                   previously using snpindex) for tolerance to SNPs

  --split-output=STRING          Basename for multiple-file output, separately for nomapping,

                                   uniq, mult, (and chimera, if --chimera-margin is selected)

  --failed-input=STRING          Print completely failed alignments as input FASTA or FASTQ format

                                   to the given file.  If the

--split-output flag is also given, this file

                                   is generated in addition to the output in the .nomapping file.

  --append-output                When --split-output or --failedinput is given, this flag will append output to the existing files.  Otherwise, the default is to create new files.

  --output-buffer-size=INT       Buffer size, in queries, for output thread (default 1000).  When the number

                                   of results to be printed exceeds this size, worker threads wait

                                   until the backlog is cleared

  --translation-code=INT         Genetic code used for translating codons to amino acids and computing CDS

                                   Integer value (default=1) corresponds to an available code at http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi

  --alt-start-codons             Also, use the alternate initiation codons shown in the above Web site

                                   By default, without this option,

only ATG is considered an initiation codon

  -F, --fulllength               Assume full-length protein, starting with Met

  -a, --cdsstart=INT             Translate codons from given

nucleotide (1-based)

  -T, --truncate                 Truncate alignment around full-length

protein, Met to Stop

                                 Implies -F flag.

  -Y, --tolerant                 Translates cDNA with corrections for

frameshifts

 

Options for GFF3 output

  --gff3-add-separators=INT      Whether to add a ### separator after each query sequence

                                   Values: 0 (no), 1 (yes, default)

  --gff3-swap-phase=INT          Whether to swap phase (0 => 0, 1 => 2, 2 => 1) in gff3_gene format

                                   Needed by some analysis programs, but deviates from GFF3 specification

                                   Values: 0 (no, default), 1 (yes)

  --gff3-fasta-annotation=INT    Whether to include annotation from the FASTA header into the GFF3 output

                                   Values: 0 (default): Do not include

                                           1: Wrap all annotation as Annot="<header>"

                                           2: Include key=value pairs, replacing brackets with quotation marks and replacing spaces between key=value pairs with semicolons

  --gff3-cds=STRING              Whether to use cDNA or genomic translation for the CDS coordinates Values: cdna (default), genomic Options for SAM output

  --no-sam-headers               Do not print headers beginning with '@'

  --sam-use-0M                   Insert 0M in CIGAR between adjacent insertions and deletions

                                   Required by Picard, but can cause errors in other tools

  --sam-extended-cigar           Use extended CIGAR format (using X and = symbols instead of M,

                                   to indicate matches and mismatches, respectively

  --force-xs-dir                 For RNA-Seq alignments, disallows XS:A:? when the sense direction

                                   is unclear, and replaces this value arbitrarily with XS:A:+.

                                   May be useful for some programs, such as Cufflinks, that cannot

                                   handle XS:A:?.  However, if you use this flag, the reported value

                                   of XS:A:+ in these cases will not be meaningful.

  --md-lowercase-snp             In MD string, when known SNPs are given by the -v flag,

                                   prints difference nucleotides as lower-case when they,

                                   differ from reference but match a known alternate allele

  --action-if-cigar-error        Action to take if there is a disagreement between CIGAR length and sequence length

                                   Allowed values: ignore, warning (default), noprint, abort

                                   Note that the noprint option does not print the CIGAR string at all if there is an error, so it may break a SAM parser

  --read-group-id=STRING         Value to put into read-group id (RG-ID) field

  --read-group-name=STRING       Value to put into read-group name (RG-SM) field

  --read-group-library=STRING    Value to put into read-group library (RG-LB) field

  --read-group-platform=STRING   Value to put into read-group library

(RG-PL) field Options for quality scores

  --quality-protocol=STRING      Protocol for input quality scores. Allowed values:

                                   illumina (ASCII 64-126) (equivalent to -J 64 -j -31)

                                   sanger   (ASCII 33-126) (equivalent to -J 33 -j 0)

                                 Default is sanger (no quality print shift)

                                 SAM output files should have quality scores in sanger protocol

 

                                 Or you can specify the print shift with this flag:

  -j, --quality-print-shift=INT  Shift FASTQ quality scores by this amount in output

                                   (default is 0 for sanger protocol; to change Illumina input

                                   to Sanger output, select -31) External map file options

  -M, --mapdir=directory         Map directory

  -m, --map=iitfile              Map file.  If argument is '?' (with the quotes),

                                   this lists available map files.

  -e, --mapexons                 Map each exon separately

  -b, --mapboth                  Report hits from both strands of genome

  -u, --flanking=INT             Show flanking hits (default 0)

  --print-comment                Show comment line for each hit

 

Alignment output options

  --nolengths                    No intron lengths in alignment

  --nomargin                     No left margin in GMAP standard output (with the -A flag)

  -I, --invertmode=INT           Mode for alignments to genomic (-) strand:

                                   0=Don't invert the cDNA (default)

                                   1=Invert cDNA and print genomic (-) strand

                                   2=Invert cDNA and print genomic (+) strand

  -i, --introngap=INT            Nucleotides to show on each end of intron (default 3)

  -l, --wraplength=INT           Wrap length for alignment (default 50) Filtering output options

  --min-trimmed-coverage=FLOAT   Do not print alignments with trimmed coverage less

                                   this value (default=0.0, which means no filtering)

                                   Note that chimeric alignments will be output regardless

                                   of this filter

  --min-identity=FLOAT           Do not print alignments with identity less

                                   this value (default=0.0, which means no filtering)

                                   Note that chimeric alignments will be output regardless

                                   of this filter

 

Help options

  --check                        Check compiler assumptions

  --version                      Show version

  --help                         Show this help message

 

 

 

> gmap_build

-k flag not specified, so building main hash table with default 15-mers

-j flag not specified, so building regional hash tables with default 6-mers

 

gmap_build: Builds a gmap database for a genome to be used by GMAP or GSNAP.

Part of GMAP package, version 2021-08-25.

 

Usage: gmap_build [options...] -d <genome> [-c <transcriptome> -T

<transcript_fasta>] <genome_fasta_files>

 

You are free to name <genome> and <transcriptome> as you wish.  You

will use the same names when performing alignments subsequently using

GMAP or GSNAP.

 

Note: If adding a transcriptome to an existing genome, then there is

no need to specify the genome_fasta_files.  This way you can add

transcriptome information to an existing genome database.

 

Options:

    -D, --dir=STRING          Destination directory for installation

(defaults to gmapdb

                                directory specified at configure time)

    -d, --genomedb=STRING     Genome name (required)

 

    -n, --names=STRING        Substitute names for contigs, provided in a file.

        The file can have two formats:

 

        1.  A file with one column per line, with each line

            corresponding to a FASTA file, in the order given to

            gmap_build.  The chromosome name for each FASTA file will

            be replaced with the desired chromosome name in the file.

            Every chromosome in the FASTA must have a corresponding line

            in the file.  This is useful if you want to rename chromosomes

            with a systematic numbering pattern.

 

        2.  A file with two columns per line, separated by white

            space.  In each line, the original FASTA chromosome name

            should be in column 1 and the desired chromosome name

            will be in column 2.

 

            The meaning of file format 2 depends on whether

            --limit-to-names is specified.  If so, the genome build will

            be limited to those chromosomes in this file.  Otherwise,

            all chromosomes in the FASTA file will be included,

            but only those chromosomes in this file will be re-named, which

            provides an easy way to change just a few chromosome names.

 

        This file can be combined with the --sort=names option, in

        which the order of chromosomes is that given in the file.  In

        this case, every chromosome must be listed in the file, and

        for chromosome names that should not be changed, column 2 can

        be blank (or the same as column 1).  The option of a blank

        column 2 is allowed only when specifying --sort=names,

        because otherwise, the program cannot distinguish between a

        1-column and 2-column names file.

 

    -L, --limit-to-names      Determines whether to limit the genome

build to the lines listed

                              in the --names file.  You can limit a

genome build to certain

                              chromosomes with this option, plus a

--names file that either

                              renames chromosomes, or lists the same

names in both columns for

                              the desired chromosomes.

 

    -k, --kmer=INT            k-mer value for genomic index (allowed:

15 or less, default is 15)

    -q INT                    sampling interval for genomoe (allowed:

1-3, default 3)

 

    -s, --sort=STRING         Sort chromosomes using given method:

                                none - use chromosomes as found in

FASTA file(s) (default)

                                alpha - sort chromosomes

alphabetically (chr10 before chr 1)

                                numeric-alpha - chr1, chr1U, chr2,

chrM, chrU, chrX, chrY

                                chrom - chr1, chr2, chrM, chrX, chrY,

chr1U, chrU

                                names - sort chromosomes based on file

provided to --names flag

 

    -g, --gunzip              Files are gzipped, so need to gunzip

each file first

    -E, --fasta-pipe=STRING   Interpret argument as a command, instead

of a list of FASTA files

    -Q, --fastq               Files are in FASTQ format

    -R, --revcomp             Reverse complement all contigs

    -w INT                    Wait (sleep) this many seconds after

each step (default 2)

 

    -o, --circular=STRING     Circular chromosomes (either a list of

chromosomes separated

                                by a comma, or a filename containing

circular chromosomes,

                                one per line).  If you use the --names

feature, then you

                                should use the substitute name of the

chromosome, not the

                                original name, for this option.

(NOTE: This behavior is different

                                from previous versions, and starts

with version 2020-10-20.)

 

    -2, --altscaffold=STRING  File with alt scaffold info, listing

alternate scaffolds,

                                one per line, tab-delimited, with the

following fields:

                                (1) alt_scaf_acc, (2) parent_name, (3)

orientation,

                                (4) alt_scaf_start, (5) alt_scaf_stop,

(6) parent_start, (7) parent_end.

 

    -e, --nmessages=INT       Maximum number of messages (warnings,

contig reports) to report (default 50)

 

 

Options for older genome formats:

    -M, --mdflag=STRING       Use MD file from NCBI for mapping contigs to

                                chromosomal coordinates

 

    -C, --contigs-are-mapped  Find a chromosomal region in each FASTA

header line.

                                Useful for contigs that have been mapped

                                to chromosomal coordinates.  Ignored

if the --mdflag is provided.

 

 

Options for transcriptome-guided alignment:

    -c, --transcriptomedb=STRING    Transcriptome name

 

    -T, --transcripts=FILE    FASTA file containing transcripts

(required if specifying

                                --transcriptomedb)

 

    -t, --nthreads=INT        Number of threads for GMAP alignment of

transcripts to genome

                                (default 8)

 

 

 

実行方法

1,indexing

gmap_build -d indexDB -k 15 genome.fna

 

2、alignment

gmap -d indexDB  -B 5 -t 10 -f 3  transcripts.fasta > output.gff3
  • -B   Batch mode (default = 2)

  • -t    Number of worker threads
  • -f    Other format for output.  gff3_match_cdna (or 3) = GFF3 cDNA_match format

 

 

K-merサイズ(マニュアルより)

kフラグにより、ゲノムインデックスのk-merサイズを12から15の範囲で制御できる。kのデフォルト値は15だが、インデックスを構築するために4GBのRAMが必要になる。もし4GBのRAMがない場合は、-kの値を小さくするか、他のマシンを探す必要がある。以下は、様々なインデックスを構築するために必要なRAM。

12-mer: 64 MB
13-mer: 256 MB
14-mer: 1 GB
15-mer: 4 GB

ゲノムインデックスは圧縮されているため、インデックスの構築後にGMAP/GSNAPプログラムを実行するためのRAM容量は必要ない。例えば、k-merが15の場合のゲノムインデックスは、gammaptrsファイルが64MB、offsetscompファイルが約350MBとなり、4GBよりはるかに小さくなる。gmap_build の -b または --basesize パラメータで圧縮量を制御することができる。デフォルトでは、k-mer sizeの値は15、basesizeの値は12になっている。k-merサイズに別の値を指定すると、basesizeはそのk-merサイズと等しくなるようにデフォルトで設定される。もし、複数のk-mer sizeでゲノムデータベースを構築したい場合は、-kの値を変えてgmap_buildを再実行する。この場合、前回実行したものと同一のファイルのみが上書きされる。その後、GMAPまたはGSNAPの実行時に-kフラグを使用してk-merサイズを選択することができる。

 

 

引用

GMAP and GSNAP for Genomic Sequence Alignment: Enhancements to Speed, Accuracy, and Functionality
Thomas D Wu, Jens Reeder, Michael Lawrence, Gabe Becker, Matthew J Brauer 

Methods Mol Biol. 2016;1418:283-334

 

GMAP: a genomic mapping and alignment program for mRNA and EST sequences
Thomas D Wu 1, Colin K Watanabe

Bioinformatics. 2005 May 1;21(9):1859-75

 

関連


リファレンスアセンブリにアライメントした後のリードの品質を評価する best

 

高精度なシーケンシング技術を開発するためには、プラットフォーム依存のシーケンシングエラーを理解する必要がある。bestは、高品質のリファレンスアセンブリにアライメントされたリードを取り込み、リードごとのメトリクス、サマリー統計、ゲノム区間ごとの層別メトリクスを生成する。bestは従来の手法と比較して16倍高速であることを示す。bestは、シーケンスプラットフォームの精度を向上させる開発支援に有用であることに加え、ライブラリ調製やエラー補正方法など、他の実験的要因の評価・改善にも応用できる。bestは、MITライセンスのもと、Github (github.com/google/best) で公開されているオープンソースコマンドラインユーティリティである。

 

特徴(Githubより)

  • 全体およびアラインメントごとの統計情報を収集
  • indelの長さの分布
  • 異なる経験的 Q 値のしきい値における収量
  • 特定のタイプのリードのエラー分布を簡単に調べるためのリードごとのビン統計情報
  • インターバルで指定された領域(BEDファイル、ホモポリマー領域、ウィンドウなど)に対する統計値
  • 品質スコアと経験的Q値に関する統計情報
  • マルチスレッドによる高速化

 

インストール

Github

https://github.com/google/best

git clone https://github.com/google/best.git
cd best/
cargo build --release

> target/release/best -h

$ target/release/best -h

best 0.1.0

Daniel Liu, Daniel E. Cook

Bam Error Stats Tool (best): analysis of error types in aligned reads.

 

USAGE:

    best [OPTIONS] <INPUT> <REFERENCE> <STATS_PREFIX>

 

ARGS:

    <INPUT>           Input BAM file

    <REFERENCE>       Input reference FASTA file. Can be gzipped

    <STATS_PREFIX>    Prefix for output files that contain statistics

 

OPTIONS:

    -b, --bin-types <BIN_TYPES>...

            Types of bins to use for per alignment stats

 

    -h, --help

            Print help information

 

        --intervals-bed <INTERVALS_BED>...

            Use intervals from a BED file

 

        --intervals-border <INTERVALS_BORDER>...

            Use fixed-width window border regions as intervals

 

        --intervals-hp

            Use homopolymer regions as intervals

 

        --intervals-match <INTERVALS_MATCH>...

            Use regions that match any of the specified subsequences as intervals

 

        --intervals-window <INTERVALS_WINDOW>...

            Use fixed-width windows as intervals

 

        --intervals-window-pos <INTERVALS_WINDOW_POS>...

            Use fixed-width windows with positions as intervals

 

    -n, --name-column <NAME_COLUMN>

            Add column with a specific name in CSV outputs

 

        --no-per-aln-stats

            Turn off outputting per alignment stats

 

    -t, --threads <THREADS>

            Number of threads. Will be automatically determined if this is set to 0 [default: 0]

 

    -V, --version

            Print version information

 

 

実行方法

target/release/best input.bam reference.fasta output

output.summary_yield_stats.csv

作成中

 

引用

Best: A Tool for Characterizing Sequencing Errors

Daniel Liu, Anastasiya Belyaeva,  Kishwar Shafin,  Pi-Chuan Chang,  Andrew Carroll,  Daniel E. Cook

bioRxiv,Posted December 23, 2022