macでインフォマティクス

macでインフォマティクス

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

bamの分析に使うバイオインフォマティクスのツールキット goleft

 

goleftはMIT licence下で提供されているバイオインフォマティクスのツールキット。GO言語で構築されている。

 

インストール

Github

https://github.com/brentp/goleft

リリース(リンク)からosx向けバイナリーをダウンロードできる。パスの通ったディレクトリに移動しておく。

chmod u+x goleft-osx.dms
mv goleft-osx.dms /usr/local/bin/goleft
goleft #動作確認

$ goleft

goleft Version: 0.1.17

 

covstats   : coverage stats across bams by sampling

depth      : parallelize calls to samtools in user-defined windows

depthwed   : matricize output from depth to n-sites * n-samples

indexcov   : quick coverage estimate using only the bam index

indexsplit : create regions of even coverage across bams/crams

samplename : report samplename(s) from a bam's SM tag

goleftはBiocondaでも提供されている(リンク)。

 

 

ラン

 1、covstats   bamをサンプリングしてカバレッジとインサートサイズをレポートする。

$ goleft covstats -h

coverage insert_mean insert_sd insert_5th insert_95th template_mean template_sd pct_unmapped pct_bad_reads pct_duplicate pct_proper_pair read_length bam sample

Usage: goleft [--n N] [--regions REGIONS] BAMS [BAMS ...]

 

Positional arguments:

  BAMS                   bam for which to estimate coverage

 

Options:

  --n N, -n N            number of reads to sample for length [default: 1000000]

  --regions REGIONS, -r REGIONS

                         optional bed file to specify target regions

  --help, -h             display this help and exit

 

bamを指定してランする。bedファイルで特定の領域や染色体だけ解析もできる。

goleft covstats pair.bam

カバレッジ、インサートサイズとそのSDなどが出力される。

f:id:kazumaxneo:20180213214311j:plain

cutとcolumnに渡して整形表示。

 

 

 2、depth  一定のウィンドウサイズでbamのカバレッジを計算する。

$ goleft depth -h

Usage: goleft [--windowsize WINDOWSIZE] [--maxmeandepth MAXMEANDEPTH] [--ordered] [--q Q] [--chrom CHROM] [--mincov MINCOV] [--stats] [--reference REFERENCE] [--processes PROCESSES] [--bed BED] --prefix PREFIX BAM

 

Positional arguments:

  BAM                    bam for which to calculate depth

 

Options:

  --windowsize WINDOWSIZE, -w WINDOWSIZE

                         window size in which to calculate high-depth regions [default: 250]

  --maxmeandepth MAXMEANDEPTH, -m MAXMEANDEPTH

                         windows with depth > than this are high-depth. The default reports the depth of all regions.

  --ordered, -o          force output to be in same order as input even with -p.

  --q Q, -Q Q            mapping quality cutoff [default: 1]

  --chrom CHROM, -c CHROM

                         optional chromosome to limit analysis

  --mincov MINCOV        minimum depth considered callable [default: 4]

  --stats, -s            report sequence stats [GC CpG masked] for each window

  --reference REFERENCE, -r REFERENCE

                         path to reference fasta

  --processes PROCESSES, -p PROCESSES

                         number of processors to parallelize.

  --bed BED, -b BED      optional file of positions or regions to restrict depth calculations.

  --prefix PREFIX        prefix for output files depth.bed and callable.bed

  --help, -h             display this help and exit

 

goleft depth --reference ref.fa --prefix output input.bam 

prefixで指定したファイルに指定ウィンドウサイズのカバレッジが出力される。

 

 

 3、indexcov   bam.baiからのカバレッジの超高速な推定。30サンプルを30秒で解析可能。このコマンドは論文になっています(ref.1)。

$ goleft indexcov -h

Usage: goleft --directory DIRECTORY [--includegl] [--excludepatt EXCLUDEPATT] [--sex SEX] [--chrom CHROM] [--fai FAI] BAM [BAM ...]

 

Positional arguments:

  BAM                    bam(s) or crais for which to estimate coverage

 

Options:

  --directory DIRECTORY, -d DIRECTORY

                         directory for output files

  --includegl, -e        plot GL chromosomes like: GL000201.1 which are not plotted by default

  --excludepatt EXCLUDEPATT [default: ^chrEBV$|^NC|_random$|Un_|^HLA\-|_alt$|hap\d$]

  --sex SEX, -X SEX      comma delimited names of the sex chromosome(s) used to infer sex. Set to '' if no sex chromosomes are present. [default: X,Y]

  --chrom CHROM, -c CHROM

                         optional chromosome to extract depth. default is entire genome.

  --fai FAI, -f FAI      fasta index file. Required when crais are used.

  --help, -h             display this help and exit

 

複数のbamをディレクトリに準備してランする。

goleft indexcov --directory output/ bam_data/*.bam

関係ないファイルがあるとエラーになったので、bamとbaiだけ集めたほうがよいかもしれない。htmlで結果は出力される。

f:id:kazumaxneo:20180213220113j:plain

上では2サンプルだけです。 公式のほうがわかりやすい結果を載せています(リンク)。

 

他にも、bamをほぼ等しい量のデータでN個の領域に分割するindexsplitや(計算を並列化するために使う)、SM tagからsamplenameをレポートするsamplenameがある。

 

 

引用

ref.1

Indexcov: fast coverage quality control for whole-genome sequencing.

Pedersen BS, Collins RL, Talkowski ME, Quinlan AR

Gigascience. 2017 Nov 1;6(11):1-6.

 

https://github.com/brentp/goleft