macでインフォマティクス

macでインフォマティクス

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

多機能なNGS分析ツール BBtools 其の3BBMap追加コマンド

 

BBMapの追加コマンドについて紹介します。

 

BBMap Guide

https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbmap-guide/

 

callvariants.sh

bbcms.sh

callgenes.sh

 

 

インストール

bbmapを導入する。

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

バリアントコール

> callvariants.sh

$ callvariants.sh

 

Written by Brian Bushnell

Last modified October 12, 2017

 

Description:  Calls variants from sam or bam input.

In default mode, all input files are combined and treated as a single sample.

In multisample mode, each file is treated as an individual sample,

and gets its own column in the VCF file.  Unless overridden, input file

names are used as sample names.

Please read bbmap/docs/guides/CallVariantsGuide.txt for more information,

or bbmap/pipelines/variantPipeline.sh for a usage example.

 

Usage:  callvariants.sh in=<file,file,...> ref=<file> vcf=<file>

 

Input may be sorted or unsorted.

The reference should be fasta.

 

I/O parameters:

in=<file>       Input; may be one file or multiple comma-delimited files.

list=<file>     Optional text file containing one input file per line.

                Use list or in, but not both.

out=<file>      Output variant list in native format.  If the name ends

                with .vcf then it will be vcf format.

vcf=<file>      Output variant list in vcf format.

ref=<file>      Reference fasta.  Required to display ref alleles.

                Variant calling wil be more accurate with the reference.

shist=<file>    (scorehist) Output for variant score histogram.

zhist=<file>    (zygosityhist) Output for zygosity histogram.

overwrite=f     (ow) Set to false to force the program to abort rather than

                overwrite an existing file.

extended=t      Print additional variant statistics columns.

sample=         Optional comma-delimited list of sample names.

multisample=f   (multi) Set to true if there are multiple sam/bam files,

                and each should be tracked as an individual sample.

vcf0=           Optional comma-delimited list of per-sample outputs.

                Only used in multisample mode.

bgzip=f         Use bgzip for gzip compression.

 

Processing Parameters:

prefilter=f     Use a Bloom filter to exclude variants seen fewer than

                minreads times.  Doubles the runtime but greatly reduces

                memory usage.  The results are identical.

samstreamer=t   (ss) Load reads multithreaded to increase speed.

coverage=t      (cc) Calculate coverage, to better call variants.

ploidy=1        Set the organism's ploidy.

rarity=1.0      Penalize the quality of variants with allele fraction 

                lower than this.  For example, if you are interested in

                4% frequency variants, you could set both rarity and

                minallelefraction to 0.04.  This is affected by ploidy - 

                a variant with frequency indicating at least one copy

                is never penalized.

covpenalty=0.8  (lowcoveragepenalty) A lower penalty will increase the 

                scores of low-coverage variants, and is useful for 

                low-coverage datasets.

useidentity=t   Include average read identity in score calculation.

usepairing=t    Include pairing rate in score calculation.

usebias=t       Include strand bias in score calculation.

useedist=t      Include read-end distance in score calculation.

homopolymer=t   Penalize scores of substitutions matching adjacent bases.

nscan=t         Consider the distance of a variant from contig ends when 

                calculating strand bias.

32bit=f         Set to true to allow coverage tracking over depth 65535.

strandedcov=t   Tracks per-strand ref coverage to print the DP4 field.

                Requires more memory when enabled.

callsub=t       Call substitutions.

calldel=t       Call deletions.

callins=t       Call insertions.

nopassdot=f     Use . as genotype for variations failing the filter.

 

Trimming parameters:

border=5        Trim at least this many bases on both ends of reads.

qtrim=r         Quality-trim reads on this end

                   r: right, l: left, rl: both, f: don't quality-trim.

trimq=10        Quality-trim bases below this score.

 

Realignment parameters:

realign=f       Realign all reads with more than a couple mismatches.

                Decreases speed.  Recommended for aligners other than BBMap.

unclip=f        Convert clip symbols from exceeding the ends of the

                realignment zone into matches and substitutitions.

repadding=70    Pad alignment by this much on each end.  Typically,

                longer is more accurate for long indels, but greatly

                reduces speed.

rerows=602      Use this many rows maximum for realignment.  Reads longer

                than this cannot be realigned.

recols=2000     Reads may not be aligned to reference seqments longer 

                than this.  Needs to be at least read length plus

                max deletion length plus twice padding.

msa=            Select the aligner.  Options:

                   MultiStateAligner11ts:     Default.

                   MultiStateAligner9PacBio:  Use for PacBio reads, or for

                   Illumina reads mapped to PacBio/Nanopore reads.

 

Sam-filtering parameters:

minpos=         Ignore alignments not overlapping this range.

maxpos=         Ignore alignments not overlapping this range.

minreadmapq=4   Ignore alignments with lower mapq.

contigs=        Comma-delimited list of contig names to include. These 

                should have no spaces, or underscores instead of spaces.

secondary=f     Include secondary alignments.

supplimentary=f Include supplimentary alignments.

invert=f        Invert sam filters.

 

Variant-Calling Cutoffs:

minreads=2              Ignore variants seen in fewer reads.

minqualitymax=15        Ignore variants with lower max base quality.

minedistmax=20          Ignore variants with lower max distance from read ends.

minmapqmax=0            Ignore variants with lower max mapq.

minidmax=0              Ignore variants with lower max read identity.

minpairingrate=0.1      Ignore variants with lower pairing rate.

minstrandratio=0.1      Ignore variants with lower plus/minus strand ratio.

minquality=12.0         Ignore variants with lower average base quality.

minedist=10.0           Ignore variants with lower average distance from ends.

minavgmapq=0.0          Ignore variants with lower average mapq.

minallelefraction=0.1   Ignore variants with lower allele fraction.  This

                        should be adjusted for high ploidies.

minid=0                 Ignore variants with lower average read identity.

minscore=20.0           Ignore variants with lower Phred-scaled score.

clearfilters            Reset all filters to zero.  Filter flags placed after

                        the clearfilters flag will still be applied.

 

There are additionally max filters for score, quality, mapq, allelefraction,

and identity.

 

Java Parameters:

-Xmx            This will be passed to Java to set memory usage, overriding the program's automatic memory detection.

                -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.

 

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

 

 

ORF予測

>callgenes.sh

$ callgenes.sh 

 

Written by Brian Bushnell

Last modified December 19, 2018

 

Description:  Finds orfs and calls genes in unspliced prokaryotes.

This includes bacteria, archaea, viruses, and mitochondria.

Can also predict 16S, 23S, 5S, and tRNAs.

 

Usage:  callgenes.sh in=contigs.fa out=calls.gff outa=aminos.faa

 

File parameters:

in=<file>       A fasta file.

out=<file>      Output gff file.

outa=<file>     Amino acid output.

model=<file>    A pgm file or comma-delimited list.

                If unspecified a default model will be used.

compareto=      Optional reference gff file to compare with the gene calls.

 

Other parameters:

minlen=60       Don't call genes shorter than this.

trd=f           (trimreaddescription) Set to true to trim read headers after

                the first whitespace.  Necessary for IGV.

merge=f         For paired reads, merge before calling.

detranslate=f   Output canonical nucleotide sequences instead of amino acids.

recode=f        Re-encode nucleotide sequences over called genes, leaving

                non-coding regions unchanged.

 

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

 

> analyzegenes.sh

$ analyzegenes.sh

 

Written by Brian Bushnell

Last modified September 27, 2018

 

Description:  Generates a prokaryotic gene model (.pkm) for gene calling.

Input is fasta and gff files.

The .pkm file may be used by CallGenes.

 

Usage:  analyzegenes.sh in=x.fa gff=x.gff out=x.pgm

 

File parameters:

in=<file>       A fasta file or comma-delimited list of fasta files.

gff=<file>      A gff file or comma-delimited list.  This is optional;

                if present, it must match the number of fasta files.

                If absent, a fasta file 'foo.fasta' will imply the

                presence of 'foo.gff'.

out=<file>      Output pgm file.

 

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

 

 

メタゲノムエラーコレクション

> bbcms.sh

r$ bbcms.sh 

 

Written by Brian Bushnell

Last modified July 24, 2018

 

Description:  Error corrects reads and/or filters by depth, storing

kmer counts in a count-min sketch (a Bloom filter variant).

This uses a fixed amount of memory.  The error-correction algorithm is taken

from Tadpole; with plenty of memory, the behavior is almost identical to 

Tadpole.  As the number of unique kmers in a dataset increases, the accuracy 

decreases such that it will make fewer corrections.  It is still capable

of making useful corrections far past the point where Tadpole would crash

by running out of memory, even with the prefilter flag.  But if there is

sufficient memory to use Tadpole, then Tadpole is more desirable.

 

Because accuracy declines with an increasing number of unique kmers, it can

be useful with very large datasets to run this in 2 passes, with the first 

pass for filtering only using a 2-bit filter with the flags tossjunk=t and 

ecc=f (and possibly mincount=2 and hcf=0.4), and the second pass using a 

4-bit filter for the actual error correction.

 

Usage:  bbcms.sh in=<input file> out=<output> outb=<reads failing filters>

 

Example of use in error correction:

bbcms.sh in=reads.fq out=ecc.fq bits=4 hashes=3 k=31 merge

 

Example of use in depth filtering:

bbcms.sh in=reads.fq out=high.fq outb=low.fq k=31 mincount=2 ecc=f hcf=0.4

 

Error correction and depth filtering can be done simultaneously.

 

File parameters:

in=<file>       Primary input, or read 1 input.

in2=<file>      Read 2 input if reads are in two files.

out=<file>      Primary read output.

out2=<file>     Read 2 output if reads are in two files.

outb=<file>     (outbad/outlow) Output for reads failing mincount.

outb2=<file>    (outbad2/outlow2) Read 2 output if reads are in two files.

extra=<file>    Additional comma-delimited files for generating kmer counts.

ref=<file>      If ref is set, then only files in the ref list will be used

                for kmer counts, and the input files will NOT be used for

                counts; they will just be filtered or corrected.

overwrite=t     (ow) Set to false to force the program to abort rather than

                overwrite an existing file.

 

Hashing parameters:

k=31            Kmer length, currently 1-31.

hashes=3        Number of hashes per kmer.  Higher generally reduces 

                false positives at the expense of speed; rapidly

                diminishing returns above 4.

minprob=0.5     Ignore kmers with probability of being correct below this.

memmult=1.0     Fraction of free memory to use for Bloom filter.  1.0 should

                generally work; if the program crashes with an out of memory

                error, set this lower.  You may be able to increase accuracy

                by setting it slightly higher.

cells=          Option to set the number of cells manually.  By default this

                will be autoset to use all available memory.  The only reason

                to set this is to ensure deterministic output.

seed=0          This will change the hash function used.

 

Depth filtering parameters:

mincount=0      If positive, reads with kmer counts below mincount will

                be discarded (sent to outb).

hcf=1.0         (highcountfraction) Fraction of kmers that must be at least

                mincount to pass.

requireboth=t   Require both reads in a pair to pass in order to go to out.

                When true, if either read has a count below mincount, both

                reads in the pair will go to outb.  When false, reads will

                only go to outb if both fail.

tossjunk=f      Remove reads or pairs with outermost kmer depth below 2.

(Suggested params for huge metagenomes: mincount=2 hcf=0.4 tossjunk=t)

 

Error correction parameters:

ecc=t           Perform error correction.

bits=           Bits used to store kmer counts; max count is 2^bits-1.

                Supports 2, 4, 8, 16, or 32.  16 is best for high-depth data;

                2 or 4 are for huge, low-depth metagenomes that saturate the 

                bloom filter otherwise.  Generally 4 bits is recommended for 

                error-correction and 2 bits is recommended for filtering only.

ecco=f          Error-correct paired reads by overlap prior to kmer-counting.

merge=t         Merge paired reads by overlap prior to kmer-counting, and 

                again prior to correction.  Output will still be unmerged.

smooth=3        Remove spikes from kmer counts due to hash collisions.

                The number is the max width of peaks to be smoothed; range is

                0-3 (3 is most aggressive; 0 disables smoothing).

                This also affects tossjunk.

                

 

Java Parameters:

-Xmx            This will be passed to Java to set memory usage, overriding the program's automatic memory detection.

                -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.

 

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

 

 

 

実行方法

1、callvariants.sh

bam/samからバリアントコールを行う。

1-1、リファレンスにペアエンドfastqをmappingしてbamかsamを得る。

bbmap.sh ref=ref.fa in1=pair_1.fq in2=pair_2.fq out=output.sam nodisk
samtools sort -@ 12 -O BAM output.sam > sort.bam && samtools index -@ 12 sort.bam

indexディレクトリが作られ、リードがリファレンスにアライメントされる。defaultだとマシンの全コアが計算に使用される(明示するなら"t=")。それからsortコマンドでcoordinate sortedされたbamを得る。 

 

1-2、入力のbam/samとリファレンスを指定してcallvariantsを実行する。ここでは倍数性(ploidy)2で実行。

callvariants.sh ploidy=2 in=sort.bam ref=ref.fa vcf=output.vcf
  • callsub=t        Call substitutions
  • calldel=t         Call deletions
  • callins=t         Call insertions
  • ploidy=1         Set the organism's ploidy
  • in=<file>         Input; may be one file or multiple comma-delimited files
  • vcf=<file>       Output variant list in vcf format

バクテリアゲノムだと数秒でランは終わる。 vcf=を指定時の出力はVCF v4.2に従う。

 

2、callgenes.sh

contigやゲノム配列からprokaryotesの遺伝子配列を予測する。ここでは予測した遺伝子のアミノ酸配列も出力している。

callgenes.sh in=ref.fa out=output outa=amnoacids.fa

NGSのリードから実行することもできる。

callgenes.sh in=pair_1.fq in2=pair_2.fq out=output outa=amnoacids.fa

リードそれぞれから予測するため、出力は相応に大きくなることに注意。

 

 

3、bbcms.sh 

メタゲノムシーケンシングリードのエラーコレクション。

入力のペアエンドfastqと出力のペアエンドfastqをそれぞれ指定する。

bbcms.sh in=pair_1.fq in2=pair_2.fq \
out=error_corrected_1.fq.gz out2=error_corrected_2.fq.gz \
outb=failing_mincount_1.fq.gz outb2=failing_mincount_2.fq.gz
  • in=<file>        Primary input, or read 1 input.
  • in2=<file>      Read 2 input if reads are in two files.
  • out=<file>      Primary read output.
  • out2=<file>    Read 2 output if reads are in two files.
  • outb=<file>    Output for reads failing mincount.
  • outb2=<file>  Read 2 output if reads are in two files.

 

参考

Question: bbmap callvariants - how to add sample names, how to get 0/0 alleles back?

http://callvariants.sh

 

関連