macでインフォマティクス

macでインフォマティクス

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

高速なbam処理ツール biobambam2

 

2020 4/17 help更新、インストール追記

 

<Biobambam論文(*1)より>

 SAM(Sequence Alignment / Matching)およびBAM(Binary Alignment / Matching)ファイルフォーマットは、ハイスループットシーケンシングおよび得られたデータの参照ゲノムへのアライメントによって得られたシーケンスデータを保存するための標準フォーマットとなっている。どちらの形式もSAMtoolsパッケージの一部として導入された(論文より cf. [1]を参照)。 SAMは人間が読めるテキスト形式だが、BAMはよりコンパクトで圧縮されたバイナリ形式である。フォーマットの現在の仕様は、[ref.2]で利用できる。これらのファイルは、含まれているデータとリファレンスとの間のバリエーションの検出、シーケンシングのクオリティ管理、長期保存などの多くのアプリケーションに使用できる。 SSAHA [ref.3]、BWA [ref.4,5]、Bowtie [ref.6,7]、SOAP [ref.8,9]、SMALT [ref.10]などのリファレンス配列へのシークエンシングリードのアライメントのための多くのプログラムが作成されており、公開されたアライナは、SAMまたはBAM出力を生成することができ、またはその出力をSAMまたはBAMに変換するためのスクリプトを備える。短い線状DNA鋳型は、両端から配列決定される。これにより、各テンプレートに対して一対のシーケンスが行われる。ペアの両端には、結果のデータファイル内で同じリード名がアサインされ、その結果、両端がゲノム内の特定の予想距離内にある可能性が高いという情報を提供する。この情報は、得られたショートリード配列をリファレンスに正しくアライメントさせるか、フラグメントを新しいドラフトリファレンスにアセンブリするのを助ける。シーケンサーから得られたデータでは、ペアは通常、何らかの形で照合される。リードの順序はリードとリファンレス間のバリエーションを呼び出すのに最も適したリファレンスに揃えられる。したがって、多くのSAMおよびBAMファイルがこの順序で処理される。ただし、各ペアからの完全な情報が必要なアプリケーションもある。これには、BAMファイルのFastQ形式への変換や、バリアントの代替検出法である(論文の例:[11]参照)デノボアセンブリと重複リードのマーキングが含まれる。(一部略)重複マーキングを適用するためには、同じ座標にマッピングされたペアエンドリードのクラスタを検出する必要があるので、照合後の順序を可能な限り座標のソート順に近づけることがさらに望ましい。このホワイトペーパーでは、効率的なリード名照合アルゴリズムを使用してBAMファイルを処理するツールセットであるbiobambamを紹介する。

 含まれているツールbamtofastqとbammarkduplicates2は、広く使用されているPicardスイート([ref.12]参照)の同等のツールよりも、実行時間とメモリ要件の面でより効率的である。 bamUpdilは、bamUtilの重複検出ツール([ref.13]参照)と同様のパフォーマンスを提供するが、より大規模で複雑なデータセットを簡単に処理できる。 bamtofastqはbamUtilのbam2fastqコンポーネント紹介)より高速である。これらのフロントエンドプログラムに加えて、照合アルゴリズムは、libmausプロジェクトの他のプロジェクトで使用するためのAPIとして直接アクセスできるようにする。

  

インストール

cent osでテストした。

本体 Github

https://github.com/gt1/biobambam2

ビルド済みのバイナリをダウンロードして解凍し、bin/にパスを通す。Darwin系列のosでもビルドできる。boostとlibmaus2(Github)が必要。

 追記

#bioconda (link)
conda install -c bioconda -y biobambam

 

ラン

bamsormadup: parallel sorting and duplicate marking

s# bamsormadup  -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

level=<[-1]>                                          : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

templevel=<[1]>                                       : compression setting for temporary files (see level for options)

threads=<[12]>                                        : number of threads

tmpfile=<[bamsormadup_5b3f8fb1a727_17079_1587178270]> : prefix for temporary files, default: create files in current directory

inputformat=<[bam]>                                   : input format (sam,bam)

M=<filename>                                          : metrics file, stderr if unset

seqchksumhash=<[crc32prod]>                           : seqchksum digest function: crc32prod,crc32,md5,sha1,sha224,sha256,sha384,sha512,crc32prime32,crc32prime64,md5prime64,sha1prime64,sha224prime64,sha256prime64,sha384prime64,sha512prime64,crc32prime96,md5prime96,sha1prime96,sha224prime96,sha256prime96,sha384prime96,sha512prime96,crc32prime128,md5prime128,sha1prime128,sha224prime128,sha256prime128,sha384prime128,sha512prime128,crc32prime160,md5prime160,sha1prime160,sha224prime160,sha256prime160,sha384prime160,sha512prime160,crc32prime192,md5prime192,sha1prime192,sha224prime192,sha256prime192,sha384prime192,sha512prime192,crc32prime224,md5prime224,sha1prime224,sha224prime224,sha256prime224,sha384prime224,sha512prime224,crc32prime256,md5prime256,sha1prime256,sha224prime256,sha256prime256,sha384prime256,sha512prime256,null,sha512primesums,sha512primesums512,sha3_256,murmur3,murmur3primesums128

digest=<[md5]>                                        : hash digest computed for output stream (md5, sha512)

digestfilename=<>                                   : name of file for storing hash digest computed for output stream (not stored by default)

indexfilename=<>                                    : name of file for storing BAM index (not stored by default)

SO=<coordinate|queryname>                             : output sort order (coordinate by default)

outputformat=<[bam]>                                  : output format (sam,bam,cram)

reference=<[]>                                        : reference FastA for writing cram

optminpixeldif=<[100]>                                : pixel difference threshold for optical duplicates (default: 100)

rcsupport=<[0]>                                       : add rc tag (unclipped coordinate, default: 0)

バージョン2からの新ツール。duplicateにタグをつけ、coordinateソートも行って出力する。

bamsormadup < input.bam level=1 threads=8 inputformat=bam SO=coordinate outputformat=bam indexfilename=output.bam.bai > output.bam 

処理の大半が並列化されており、I/Oがリミットになりうるため、I/Oの高速なSSDなどのストレージ使用が推奨されている。threadsフラグがなければ利用できる全ての論理コアを使う。

 

bamcollate2: reads BAM and writes BAM reordered such that alignment or collated by query name

# bamcollate2 -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

collate=<[1]>                                   : collate pairs

reset=<>                                        : reset alignments and header like bamreset (for collate=0,1 or 3 only, default enabled for 3)

level=<[-1]>                                    : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

filename=<[stdin]>                              : input filename (default: read file from standard input)

inputformat=<[bam]>                             : input format: cram, bam or sam

reference=<>                                  : name of reference FastA in case of inputformat=cram

ranges=<>                                     : input ranges (bam/cram input only, collate<2 only, default: read complete file)

exclude=<[SECONDARY,SUPPLEMENTARY]>             : exclude alignments matching any of the given flags

disablevalidation=<[0]>                         : disable validation of input data

colhlog=<[18]>                                  : base 2 logarithm of hash table size used for collation

colsbs=<[134217728]>                            : size of hash table overflow list in bytes

T=<[bamcollate2_5b3f8fb1a727_17081_1587178319]> : temporary file name

md5=<[0]>                                       : create md5 check sum (default: 0)

md5filename=<filename>                          : file name for md5 check sum (default: extend output file name)

index=<[0]>                                     : create BAM index (default: 0)

indexfilename=<filename>                        : file name for BAM index file (default: extend output file name)

readgroups=[<>]                                 : read group filter (default: keep all)

mapqthres=<[-1]>                                : mapping quality threshold (collate=1 only, default: keep all)

classes=[F,F2,O,O2,S]                           : class filter (collate=1 only, default: keep all)

resetheadertext=[<>]                            : replacement SAM header text file for reset=1 (default: filter header in source BAM file)

resetaux=<[1]>                                  : reset auxiliary fields (collate=0,1 only with reset=1)

auxfilter=[<>]                                  : comma separated list of aux tags to keep if reset=1 and resetaux=0 (default: keep all)

outputformat=<[bam]>                            : output format (bam,cram,sam)

O=<[stdout]>                                    : output filename (standard output if unset)

outputthreads=<[1]>                             : output helper threads (for outputformat=bam only, default: 1)

inputbuffersize=<[65536]>                       : input buffer size

replacereadgroupnames=<[]>                      : name of file containing list of read group identifier replacements (default: no replacements)

 

Alignment flags: PAIRED,PROPER_PAIR,UNMAP,MUNMAP,REVERSE,MREVERSE,READ1,READ2,SECONDARY,QCFAIL,DUP,SUPPLEMENTARY

ペアエンドを照合し、secondary hitは除く。samtools collateと同様の機能(samtools1.8 manual)。

bamcollate2 <input.bam collate=1 level=1 exclude=SECONDARY > out.bam 

 level=0だと非圧縮になる。

 

bammarkduplicates: reads BAM and writes BAM with duplicate alignments marked using the BAM flags field

# bammarkduplicates -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

I=<filename>                  : input file, stdin if unset

O=<filename>                  : output file, stdout if unset

M=<filename>                  : metrics file, stderr if unset

D=<filename>                  : duplicates output file if rmdup=1

tmpfile=<filename>            : prefix for temporary files, default: create files in current directory

level=<[-1]>                  : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

markthreads=<[1]>             : number of helper threads

verbose=<[1]>                 : print progress report (default: 1)

mod=<[1048576]>               : print progress for each mod'th record/alignment

rewritebam=<[0]>              : compression of temporary alignment file when input is via stdin (0=snappy,1=gzip/bam,2=copy)

rewritebamlevel=<[-1]>        : compression setting for rewritten input file if rewritebam=1 (-1=zlib default,0=uncompressed,1=fast,9=best)

rmdup=<[0]>                   : remove duplicates (default: 0)

md5=<[0]>                     : create md5 check sum (default: 0)

md5filename=<filename>        : file name for md5 check sum (default: extend output file name)

index=<[0]>                   : create BAM index (default: 0)

indexfilename=<filename>      : file name for BAM index file (default: extend output file name)

tag=<[a-zA-Z][a-zA-Z0-9]>     : aux field id for tag string extraction

nucltag=<[a-zA-Z][a-zA-Z0-9]> : aux field id for nucleotide tag extraction

colhashbits=<[20]>            : log_2 of size of hash table used for collation

collistsize=<[33554432]>      : output list size for collation

fragbufsize=<[50331648]>      : size of each fragment/pair file buffer in bytes

D=<filename>                  : duplicates output file if rmdup=1

dupmd5=<[0]>                  : create md5 check sum for duplicates output file (default: 0)

dupmd5filename=<filename>     : file name for md5 check sum of dup file (default: extend duplicates output file name)

dupindex=<[0]>                : create BAM index for duplicates file (default: 0)

dupindexfilename=<filename>   : file name for BAM index file for duplicates file (default: extend duplicates output file name)

optminpixeldif=<[100]>        : pixel difference threshold for optical duplicates (default: 100)

dupilcateを除く。

./bammarkduplicates <input.bam rmdup=1 markthreads=8 level=1 > out.bam 

 

 

bammaskflags: reads BAM and writes BAM while masking (removing) bits from the flags column

# bammaskflags -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

level=<[-1]>             : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

maskneg=<[0]>            : flag masking bitmask (bits set are removed)

md5=<[0]>                : create md5 check sum (default: 0)

md5filename=<filename>   : file name for md5 check sum (default: extend output file name)

index=<[0]>              : create BAM index (default: 0)

indexfilename=<filename> : file name for BAM index file (default: extend output file name)

tmpfile=<filename>       : prefix for temporary files, default: create files in current directory

auxexists=[]             : only mask records containing any of the given aux field tags

 

該当するflag情報(sam/bamの2列目)をもつリードを除く。

./bammaskflags < input.bam level=1 maskneg=0 > out.bam 

 

bamrecompress: reads BAM and writes BAM with a defined compression setting. This tool is capable of multi-threading.

# bamrecompress -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

level=<[-1]>             : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

verbose=<[1]>            : print progress report

numthreads=<[1]>         : number of recoding threads

md5=<[0]>                : create md5 check sum (default: 0)

md5filename=<filename>   : file name for md5 check sum (default: extend output file name)

index=<[0]>              : create BAM index (default: 0)

indexfilename=<filename> : file name for BAM index file (default: extend output file name)

tmpfile=<filename>       : prefix for temporary files, default: create files in current directory

bamを再圧縮する。

bamrecompress <input.bam level=9 numthreads=8 > output.bam

 

 

bamsort: reads BAM and writes BAM resorted by coordinates or query name

# bamsort -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

level=<[-1]>                  : compression settings for output bam file (-1=zlib default,0=uncompressed,1=fast,9=best)

SO=<[coordinate]>             : sorting order (coordinate, queryname, hash, tag or queryname_HI)

verbose=<[1]>                 : print progress report

blockmb=<[1024]>              : size of internal memory buffer used for sorting in MiB

disablevalidation=<[0]>       : disable input validation (default is 0)

tmpfile=<filename>            : prefix for temporary files, default: create files in current directory

md5=<[0]>                     : create md5 check sum (default: 0)

md5filename=<filename>        : file name for md5 check sum

index=<[0]>                   : create BAM index (default: 0)

indexfilename=<filename>      : file name for BAM index file

inputformat=<[bam]>           : input format (bam,cram,maussam,sam,sbam)

outputformat=<[bam]>          : output format (bam,cram,sam)

I=<[stdin]>                   : input filename (standard input if unset)

inputthreads=<[1]>            : input helper threads (for inputformat=bam only, default: 1)

reference=<>                  : reference FastA (.fai file required, for cram i/o only)

range=<>                      : coordinate range to be processed (for coordinate sorted indexed BAM input only)

outputthreads=<[1]>           : output helper threads (for outputformat=bam only, default: 1)

O=<[stdout]>                  : output filename (standard output if unset)

fixmates=<[0]>                : fix mate information (for name collated input only, disabled by default)

calmdnm=<[0]>                 : calculate MD and NM aux fields (for coordinate sorted output only)

calmdnmreference=<>         : reference for calculating MD and NM aux fields (calmdnm=1 only)

calmdnmrecompindetonly=<[0]>  : only recalculate MD and NM in the presence of indeterminate bases (calmdnm=1 only)

calmdnmwarnchange=<[0]>       : warn when changing existing MD/NM fields (calmdnm=1 only)

adddupmarksupport=<[0]>       : add info for streaming duplicate marking (for name collated input only, ignored for fixmate=0, disabled by default)

tag=<[a-zA-Z][a-zA-Z0-9]>     : aux field id for tag string extraction (adddupmarksupport=1 only)

nucltag=<[a-zA-Z][a-zA-Z0-9]> : aux field id for nucleotide tag extraction (adddupmarksupport=1 only)

markduplicates=<[0]>          : mark duplicates (only when input name collated and output coordinate sorted, disabled by default)

rmdup=<[0]>                   : remove duplicates (only when input name collated and output coordinate sorted, disabled by default)

streaming=<[1]>               : do not open input files multiple times when set

sorttag=<>                  : tag used by SO=tag (no default)

sortthreads=<[1]>             : threads used for sorting (default: 1)

 

ポジションソートしてbam出力する。

./bamsort SO=coordinate < input.bam inputthreads=8 level=1 outputformat=bam > out.bam 

 

bamtofastq: reads BAM and writes FastQ; output can be collated or uncollated by query name

# bamtofastq -h

This is biobambam2 version 2.0.87.

biobambam2 is distributed under version 3 of the GNU General Public License.

 

Key=Value pairs:

 

F=<[stdout]>                                   : matched pairs first mates

F2=<[stdout]>                                  : matched pairs second mates

S=<[stdout]>                                   : single end

O=<[stdout]>                                   : unmatched pairs first mates

O2=<[stdout]>                                  : unmatched pairs second mates

collate=<[1]>                                  : collate pairs

combs=<[0]>                                    : print some counts after collation based processing

filename=<[stdin]>                             : input filename (default: read file from standard input)

inputformat=<[bam]>                            : input format: cram, bam or sam

reference=<>                                 : name of reference FastA in case of inputformat=cram

ranges=<>                                    : input ranges (bam and cram input only, default: read complete file)

exclude=<[SECONDARY,SUPPLEMENTARY]>            : exclude alignments matching any of the given flags

disablevalidation=<[0]>                        : disable validation of input data

colhlog=<[18]>                                 : base 2 logarithm of hash table size used for collation

colsbs=<[33554432]>                            : size of hash table overflow list in bytes

T=<[bamtofastq_5b3f8fb1a727_17091_1587178446]> : temporary file name

gz=<[0]>                                       : compress output streams in gzip format (default: 0)

level=<[-1]>                                   : compression setting if gz=1 (-1=zlib default,0=uncompressed,1=fast,9=best)

fasta=<[0]>                                    : output FastA instead of FastQ

inputbuffersize=<[-1]>                         : size of input buffer

outputperreadgroup=<[0]>                       : split output per read group (for collate=1 only)

outputdir=<>                                   : directory for output if outputperreadgroup=1 (default: current directory)

outputperreadgroupsuffixF=<[_1.fq]>            : suffix for F category when outputperreadgroup=1

outputperreadgroupsuffixF2=<[_2.fq]>           : suffix for F2 category when outputperreadgroup=1

outputperreadgroupsuffixO=<[_o1.fq]>           : suffix for O category when outputperreadgroup=1

outputperreadgroupsuffixO2=<[_o2.fq]>          : suffix for O2 category when outputperreadgroup=1

outputperreadgroupsuffixS=<[_s.fq]>            : suffix for S category when outputperreadgroup=1

tryoq=<[0]>                                    : use OQ field instead of quality field if present (collate={0,1} only)

split=<[0]>                                    : split named output files into chunks of this amount of reads (0: do not split)

splitprefix=<[bamtofastq_split]>               : file name prefix if collate=0 and split>0

tags=<>                                      : list of aux tags to be copied (default: do not copy any aux fields)

outputperreadgrouprgsm=<[0]>                   : add read group field SM ahead of read group id when outputperreadgroup=1 (for collate=1 only)

outputperreadgroupprefix=<>                  : prefix added in front of file names if outputperreadgroup=1 (for collate=1 only)

 

Alignment flags: PAIRED,PROPER_PAIR,UNMAP,MUNMAP,REVERSE,MREVERSE,READ1,READ2,SECONDARY,QCFAIL,DUP,SUPPLEMENTARY

bamからペアエンドfastqを出力。

bamtofastq < out.bam F=1.fq F2=2.fq

 

 

引用

*1

biobambam: tools for read pair collation based algorithms on BAM files

German Tischler, Steven Leonard

Source Code Biol Med. 2014; 9: 13.