macでインフォマティクス

macでインフォマティクス

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

複数のリファレンスに同時にマッピングしてリードを分ける BBSplit

 

SEQanswersより。一部改変(リンク

BBSplitは、BBMapを用いて、複数のリファレンスに同時にマッピングすることで、リードをビン分け(binning)します。リードは、最もよくマップされるリファレンスのビンに分けて書き出されます。また、曖昧さ回避のオプションもあり、複数のリファレンスにマッピングされたリードは、それらのすべて、なし、いずれか、またはそれぞれの特別な「ambiguous」ファイルに入れられるようにできます。ペアエンドのリードは同期されて保存されます。

使用例として、例えば、大腸菌サルモネラ菌に汚染された何かのライブラリを持っていた場合、次のようにすることができます。

bbsplit.sh in1=reads1.fq in2=reads2.fq ref=ecoli.fa,salmonella.fa basename=out_%.fq outu1=clean1.fq outu2=clean2.fqのようにする。

これにより、3つの出力ファイルが生成されます。

パラメータに関する詳しい情報は、bbsplit.shを引数なしで実行するか、/bbmap/docs/readme.txtを参照して下さい。

OSがbashシェルスクリプトを処理できない場合は、「bbsplit.sh」を”java -Xmx29g -cp /path/to/current align2.BBSplitter”に置き換えてください。ここで/path/to/currentは「現在の」ディレクトリ(bbmapのサブディレクトリ)の場所で、-Xmx29gは使用するメモリ量(つまりこれは32GBコンピュータに対するコマンドラインとなる)を指定します。これは物理メモリの約85%に設定する必要があります。

BBSplitはマッピングにBBMapを使用し、非常に高速かつ高感度です。そのため、BBMapでサポートされているすべてのフラグと機能はBBSplitで使用できます(sam出力は別として)。

BBSplitはこちらから入手可能です。
https://sourceforge.net/projects/bbmap/

 

BOL: BBSplit: Read Binning Tool for Metagenomes and Contaminated Librariesより

BBSplitは内部でBBMapを使用して、一度に複数のゲノムにリードをマッピングし、どのゲノムに最もマッチするかを判断する。これは、通常のマッピングとは異なる。あるゲノム(例えばヒト)のどこかに正確なリピートがある場合、そこにマッピングされたリードは曖昧にしかマッピングされない。しかし、マウスかヒトかを判断する場合、リードがヒトの中で曖昧にマッピングされるかどうかは問題ではなく、ヒトとマウスの間で曖昧にマッピングされるかどうかだけが問題となる。BBSplitはこの追加的なambiguity情報を追跡し、ambig2フラグに基づいてその使い方を決定する。

 

ヒトとマウスのリードを分けたり、汚染除去などに使えるコマンドです。簡単に紹介します。

インストール

condaを使ってBBToolsを導入すれば使える。

Github

mamba install -c bioconda bbmap -y

$ bbsplit.sh 

BBSplit

Written by Brian Bushnell, from Dec. 2010 - present

Last modified June 11, 2018

 

Description:  Maps reads to multiple references simultaneously.

Outputs reads to a file for the reference they best match, with multiple options for dealing with ambiguous mappings.

 

To index:     bbsplit.sh build=<1> ref_x=<reference fasta> ref_y=<another reference fasta>

To map:       bbsplit.sh build=<1> in=<reads> out_x=<output file> out_y=<another output file>

 

To be concise, and do everything in one command:

bbsplit.sh ref=x.fa,y.fa in=reads.fq basename=o%.fq

 

that is equivalent to

bbsplit.sh build=1 in=reads.fq ref_x=x.fa ref_y=y.fa out_x=ox.fq out_y=oy.fq

 

By default paired reads will yield interleaved output, but you can use the # symbol to produce twin output files.

For example, basename=o%_#.fq will produce ox_1.fq, ox_2.fq, oy_1.fq, and oy_2.fq.

 

     

Indexing Parameters (required when building the index):

ref=<file,file>     A list of references, or directories containing fasta files.

ref_<name>=<ref.fa> Alternate, longer way to specify references. e.g., ref_ecoli=ecoli.fa

                    These can also be comma-delimited lists of files; e.g., ref_a=a1.fa,a2.fa,a3.fa

build=<1>           If multiple references are indexed in the same directory, each needs a unique build ID.

path=<.>            Specify the location to write the index, if you don't want it in the current working directory.

 

Input Parameters:

build=<1>           Designate index to use.  Corresponds to the number specified when building the index.

in=<reads.fq>       Primary reads input; required parameter.

in2=<reads2.fq>     For paired reads in two files.

qin=<auto>          Set to 33 or 64 to specify input quality value ASCII offset.

interleaved=<auto>  True forces paired/interleaved input; false forces single-ended mapping.

                    If not specified, interleaved status will be autodetected from read names.

 

Mapping Parameters:

maxindel=<20>       Don't look for indels longer than this.  Lower is faster.  Set to >=100k for RNA-seq.

minratio=<0.56>     Fraction of max alignment score required to keep a site.  Higher is faster.

minhits=<1>         Minimum number of seed hits required for candidate sites.  Higher is faster.

ambiguous=<best>    Set behavior on ambiguously-mapped reads (with multiple top-scoring mapping locations).

                       best   (use the first best site)

                       toss   (consider unmapped)

                       random   (select one top-scoring site randomly)

                       all   (retain all top-scoring sites.  Does not work yet with SAM output)

ambiguous2=<best>   Set behavior only for reads that map ambiguously to multiple different references.

                    Normal 'ambiguous=' controls behavior on all ambiguous reads;

                    Ambiguous2 excludes reads that map ambiguously within a single reference.

                       best   (use the first best site)

                       toss   (consider unmapped)

                       all   (write a copy to the output for each reference to which it maps)

                       split   (write a copy to the AMBIGUOUS_ output for each reference to which it maps)

qtrim=<true>        Quality-trim ends to Q5 before mapping.  Options are 'l' (left), 'r' (right), and 'lr' (both).

untrim=<true>       Undo trimming after mapping.  Untrimmed bases will be soft-clipped in cigar strings.

 

Output Parameters:

out_<name>=<file>   Output reads that map to the reference <name> to <file>.

basename=prefix%suffix     Equivalent to multiple out_%=prefix%suffix expressions, in which each % is replaced by the name of a reference file.

bs=<file>           Write a shell script to 'file' that will turn the sam output into a sorted, indexed bam file.

scafstats=<file>    Write statistics on how many reads mapped to which scaffold to this file.

refstats=<file>     Write statistics on how many reads were assigned to which reference to this file.

                    Unmapped reads whose mate mapped to a reference are considered assigned and will be counted.

nzo=t               Only print lines with nonzero coverage.

 

***** Notes *****

Almost all BBMap parameters can be used; run bbmap.sh for more details.

Exceptions include the 'nodisk' flag, which BBSplit does not support.

BBSplit is recommended for fastq and fasta output, not for sam/bam output.

When the reference sequences are shorter than read length, use Seal instead of BBSplit.

 

Java Parameters:

-Xmx                This will set Java's memory usage, overriding autodetection.

                    -Xmx20g will specify 20 gigs of RAM, and -Xmx200m will specify 200 megs.

                    The max is typically 85% of physical memory.

-eoom               This flag will cause the process to exit if an

                    out-of-memory exception occurs.  Requires Java 8u92+.

-da                 Disable assertions.

 

This list is not complete.  For more information, please consult /readme.txt

Please contact Brian Bushnell at bbushnell@lbl.gov if you encounter any problems.

 

 

実行方法

ペアエンドのリードと、2つ以上のリファレンス配列(fasta形式)を指定する。

bbsplit.sh in1=reads1.fq.gz in2=reads2.fq.gz ref=ecoli.fa,salmonella.fa basename=out_%.fq.gz outu1=clean1.fq.gz outu2=clean2.fq.gz
  • ambiguous=<best>    Set behavior on ambiguously-mapped reads (with multiple top-scoring mapping locations).
                           best   (use the first best site)
                           toss   (consider unmapped)
                           random   (select one top-scoring site randomly)
                           all   (retain all top-scoring sites.  Does not work yet with SAM output)
  • ambiguous2=<best>   Set behavior only for reads that map ambiguously to multiple different references. Normal 'ambiguous=' controls behavior on all ambiguous reads; Ambiguous2 excludes reads that map ambiguously within a single reference.
                           best   (use the first best site)
                           toss   (consider unmapped)
                           all   (write a copy to the output for each reference to which it maps)
                           split   (write a copy to the AMBIGUOUS_ output for each reference to which it maps)

リファレンスの配列サイズがヒトゲノムのように大きいと、マッピング前のindexingに数分ほど時間がかかる。

 

モードオプションは、複数のリファレンスにマップされるリードと、それぞれのリードが異なるリファレンスにマップされるペアをどうするかを決定します。

ambig2=best
デフォルト。曖昧なリードは最初の最適なサイトへ行く。

ambig2=toss
曖昧なリードはマッピングされていないと見なす。

ambig2=all
マッピングされる各リファレンスのコピーを出力に書き込む。

ambig2=split
AMBIGUOUS_出力ファイルに、マッピングされた各参照のコピーを書き込む。

 

出力

それぞれのリファレンスにマッピングされたペアエンドリード(インターレース)と、どのリファレンスにもマッピングされなかったリード(clean~)がそれぞれ出力される。ペアエンドはいずれの出力でも同期されている。

f:id:kazumaxneo:20220313211937p:plain

必要なら、インターレースのペアリードデータを2つのファイルに分ける。

reformat.sh in=interlace.fq.gz out1=out_R1.fq.gz out2=out_R2.fq.gz

 

 

それぞれのリファレンスゲノムのどの染色体にどれだけリードがマッピングされたかの定量結果をプリントするにはscafstats=を付ける。

reformat.sh in=interlace.fq.gz out1=out_R1.fq.gz out2=out_R2.fq.gz scafstats=mapping.txt
  • scafstats=<file>    Write statistics on how many reads mapped to which scaffold to this file

 

2つ以上指定したリファレンスのどれにリードがマッピングされたかの定量結果をプリントするにはscafstats=を付ける。

reformat.sh in=interlace.fq.gz out1=out_R1.fq.gz out2=out_R2.fq.gz scafstats=mapping.txt
  • refstats=<file>     Write statistics on how many reads were assigned to which reference to this file. Unmapped reads whose mate mapped to a reference are considered assigned and will be counted.

 

引用

http://seqanswers.com/forums/showthread.php?t=41288

Gihtub

https://github.com/BioInfoTools/BBMap/blob/master/sh/bbsplit.sh

 

https://bioinformaticsonline.com/pages/view/35033/bbsplit-read-binning-tool-for-metagenomes-and-contaminated-libraries

 

関連