macでインフォマティクス

macでインフォマティクス

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

germlineとsomaticの変異を検出する SNVSniffer

 

 次世代シークエンシング(NGS)に基づいて、単一ヌクレオチド変異(SNV)または挿入 - または欠失(indel)突然変異を呼び出すための様々なアプローチが開発されている。しかし、それらの大部分は、特定のタイプの突然変異に捧げられている。正常細胞における生殖系列SNV、癌/腫瘍細胞における体細胞SNV、またはindelsのみである。文献では、生殖系列および体細胞SNVs / indelsの両方に対する効率的かつ統合されたコーラーは、まだ広範に調査されていない。

 SNVSnifferは、NGSデータからの生殖系列と体細胞SNVs / indelsの両方を同定する効率的かつ統合されたコーラーである。このアルゴリズムでは、ベイジアン確率モデルを使用してSNVを特定し、indelsを呼び出すため複数のungappedアライメントアプローチを調べることを提案する(一部略)。

 多様なベンチマークを使用して、包括的なパフォーマンス評価を実施した。生殖系列変異コールのために、SNVSnifferは最先端のFaSD、GATKおよびSAMtoolsと比較して優れたスピードで競争力の高い正確性を示す。体細胞変異のコールのために、このアルゴリズムは、VarScan2、SomaticSniper、JointSNVMix2、MuTectよりも高速で、同等またはそれ以上の精度を達成する。

 

 

本ツールには、germlineとsomaticの変異を再現可能なシーケンスリードのシミュレータ機能も備わっています。ユーザー定義の条件で変異を発生させて、本ツールで解析し、評価する、という流れが本ツールのみで可能になっています。

 

公式サイト 

http://snvsniffer.sourceforge.net/homepage.htm#latest

 version2マニュアル

http://snvsniffer.sourceforge.net/version2.htm#latest

 

インストール

cent OSに導入した。

依存

  • SAMtools

 

本体は上記公式から実行ファイルをダウンロードして使う。SNVSniffer 1.0はSNVのみがターゲットだが、SNVSniffer 2.0 はindelも検出可能になっている。

$ ./SNVSniffer 

 

SNVSniffer (v2.0.4): identifying germline and somatic single nucleotide and indel mutations (released on March 07, 2016)

Command:

snp            Identify single-sample SNVs/indels

somatic        Identify germline/somatic SNVs/indels from tumor-normal pairs

gsim           Simulate sample reads with germline SNVs/indels

ssim           Simulate normal and tumor sample reads with somatic SNVs/indels

eval           Evaluate predicted germline/somatic SNVs/indels against gold standard, only for evaluation purpose

 

  

 

ラン

5つのコマンドのうち、performance結果を評価するevalコマンド以外の4つを紹介する。テストデータがアップされているので、それを使ってテストランする。リファレンスは"Home / data / v1.0 / synthetic_tumors_error_0.01"にある(リンク)。

 

1、snp - Identify single-sample SNVs/indels

$ SNVSniffer snp

 

Usage: SNVSniffer snp [options] sam_header infile

* sam_header: a SAM header file associated with the BAM file

* infile: a mpileup/pileup/BAM file (can have no header)

Input:

-g <string> (reference genome file, [Required by BAM format])

-f <int> (input file format, default = 2)

    0: mpileup format generated by SAMtools

    1: pipeup format generated by MAQ

    2: BAM file

 

Output:

-o <string> (output file name, default = STDOUT)

 

Base call and coverage:

-min_cov <int> (minimum coverage, default = 3)

-max_cov <int> (maximum coverage, default = 65536)

-min_bqual <int> (minimum Phred base quality score, default = 13)

-min_mapq <int> (minimum mapping quality score, default = 0)

-seq_err_rate <float> (sequencing error rate of the data, default = 0.01)

 

SNV calling:

-exe_mode <int> (execution model for SNV calling, default = 0)

    0: very fast and lesst accurate

    1: fast and accurate

    2: slow and most accurate

-prior <int> (model used for genotype prior probabilities, default = 1)

    0: equal probability

    1: prior probabilities with no consideration of Ti/Tv ratio

    2: prior probabilties by considering Ti/Tv ratio

-snp_rate <float> (SNP mutation rate for the species, default = 0.001)

-min_allele_freq <int> (minimum allele frequency, default = 0.2 [Important])

-stringency <int> (stringency level [0, 9] for low-confidence mutations, default = 6)

-homo_freq <float> (allelic frequency threshold for homozygous genotype, default = 0.75)

-relaxed <int> (using the relaxed confidence score calcuation (favours sensitivity), default = 1)

-min_locus_dist <int> (minimum distance interval between two neighboring SNP loci, default = 1 [0 disables it])

-use_strand_dist <bool> (use strand distribution information, default = 0)

-local_range <int> (set the local search range for strand distribution, default = 1000)

-call_snps <bool> (enable the calling of SNPs, default = 1)

-call_indel <bool> (enable the calling of insertions or deletions, default = 1)

-pvalue <float> (P value threshold for variant calling, default = 0.05)

ランにはSAMのヘッダー情報が必要。ヘッダーファイルを作成してからランする。

samtools view -H input.bam > header.sam 
SNVSniffer snp -f 2 header.sam -g reference.fa input.bam > output.vcf

出力。

$ head -n 15 out.vcf

##fileformat=VCFv4.1

##source=SNVSniffer-v2.0.4

##INFO=<ID=DP,Number=1,Type=Integer,Description="Coverage">

##INFO=<ID=HC,Number=1,Type=Float,Description="High confidence">

##INFO=<ID=LC,Number=1,Type=Float,Description="Low confidence">

##INFO=<ID=VT,Number=1,Type=String,Description="Variant type: SNP, INS or DEL">

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE

chr21 5010446 . T C . PASS DP=16;VT=SNP;HC=3.657 GT 0|1

chr21 5012778 . C G . PASS DP=26;VT=SNP;HC=4.350 GT 0|1

chr21 5013197 . C T . PASS DP=27;VT=SNP;HC=3.657 GT 0|1

chr21 5015438 . C A . PASS DP=45;VT=SNP;HC=4.350 GT 0|1

chr21 5016893 . C G . PASS DP=33;VT=SNP;HC=4.350 GT 0|1

chr21 5017812 . G A . PASS DP=31;VT=SNP;LC=3.657 GT 0|1

chr21 5018805 . G A . PASS DP=32;VT=SNP;HC=3.657 GT 0|1 

 

 

2、somatic - Identify germline/somatic SNVs/indels from tumor-normal pairs

$ SNVSniffer somatic 

 

Usage: SNVSniffer somatic [options] normal_sam_header tumor_sam_header normal tumor

1. normal_sam_header: SAM header file for normal

2. tumor_sam_header: SAM header file for tumor

3. normal: a mpileup/pileup/BAM file for normal

4. tumor: a mpileup/pileup/BAM file for tumor

Input:

-g <string> (reference genome file, required for BAM format input)

-f <int> (input file format, default = 2)

    0: mpileup format generated by SAMtools

    1: pipeup format generated by MAQ

    2: BAM format

 

Output:

-o <string> (output file name, default = STDOUT)

-o_somatic <int> (output somatic mutations, default = 1)

-o_loh <int> (output loss of heterzygosity mutations, default = 1)

-o_germline <int> (output germline mutations, default = 0)

-o_unknown <int> (output unknown type mutations, default = 1)

-o_indel <int> (output indels, default = 1)

 

Base call and coverage:

-min_cov_normal <int> (minimum normal coverage, default = 3)

-min_cov_tumor <int> (minimum tumor coverage, default = 3)

-min_bqual <int> (minimum Phred base quality score, default = 13)

-min_mapq <int> (minimum mapping quality score, default = 0)

-seq_err_rate <float> (sequencing error rate of the data, default = 0.01)

 

 

SNV calling:

-exe_mode <int> (execution model for SNV calling, default = 0)

    0: very fast

    1: fast and more accurate

    2: slow and most accurate

-prior <int> (model used for genotype prior probabilities, default = 1)

    0: equal probability

    1: prior probabilities with no consideration of Ti/Tv ratio

    2: prior probabilties by considering Ti/Tv ratio

-somatic_rate <float> (somatic mutation rate, default = 0.01)

-tumor_purity <float> (estimated purity (tumor data) of tumor sample, default = 0 [0 means AUTO])

-min_allele_freq <float> (minimum allele frequency for the normal, default = 0.2)

-stringency <int> (stringency level [0, 9] for low-confidence mutations, default = 6)

-homo_freq <float> (allelic frequency threshold for homozygous genotype, default = 0.75)

-p_value <float> (P value threshold for variant calling per sample, default = 0.05 [0 means global])

-use_strand_dist <bool> (use strand distribution information, default = 0)

-local_range <int> (set the local search range for strand distribution, default = 1000)

 

 

最初にnormalとtumorのbamからheader.samを作る。

samtools view -H chr21_e15_c30_normal.realigned.bam> normal_header.sam
samtools view -H chr21_e15_c30_tumor.realigned.bam > tumor_header.sam

somaticコマンドを走らせる。tumorの純度を"-tumor_purity <float>"で指定しない時はSNVSnifferによって推定された値が使われる。

SNVSniffer -f 2 -g chr21.fa somatic normal_header.sam tumor_header.sam chr21_e15_c30_normal.realigned.bam chr21_e15_c30_tumor.realigned.bam -o out.vcf

 

 

3、gsim - Simulate sample reads with germline SNVs/indels

$ SNVSniffer gsim

 

SNVSniffer gsim: a diaploid germline mutation simulator derived from WGSIM in SAMtools

Usage: SVNSniffer gsim [options] reference.fa reads.base

reference.fa: reference file stored in FASTA/FASTQ format

reads.base: the base file name for reads and VCF-formated mutations

Options:

Output:

-z (compress all of the output using ZLIB)

Random numbers:

-p <int> (random number seed value for generating the simulated genome, default=11)

-q <int> (random number seed value for generating simulated reads from the simulated genome, default = 47)

Reads:

-c <int> (coverage of the genome, default = 0 [>0 distables option -n])

-n <int> (number of read pairs, default = 1000000)

-1 <int> (length of the first read, default = 100)

-2 <int> (length of the second read, default = 100)

-a <float> (disgard if the fraction of ambiguous bases higher than #FLOAT, default = 0.05)

Germline:

-e <float> (sequencing base error rate, default = 0.02)

-i <int> (average insert size, i.e. outer distance, between the two ends, default = 500)

-d <int> (standard deviation of insert size, default = 50)

-r <float> (rate of mutations, including substitutions and indels, default = 0.001)

-R <float> (fraction of indels, default = 0.15)

-X <float> (probability an indel is extended, default = 0.3)

-H <float> (homozygous variant ratio, default = 0.3333)

-h (haplotype mode (all reads are sequenced from a single sequence)

カバレッジ30で、SNVとindelが含まれたgermlineのfastqをシミュレートする。

SNVSniffer gsim reference.fa output -e 0.01 -c 30

ランが終わるとoutput_1.fqとoutput_2.fqができる。下のキャプチャ画像は bamに変換してIGVで表示してみたもの。

f:id:kazumaxneo:20180412223627j:plain

 

 

4、ssim - Simulate normal and tumor sample reads with somatic SNVs/indels

$ SNVSniffer ssim

 

SNPSniffer ssim: a diaploid somatic mutation simulator

Usage: SNVSniffer ssim [options] reference.fa reads.base

reference.fa: reference file stored in FASTA/FASTQ format

reads.base: the base file name for reads and VCF-formated mutations

Options:

Output:

-z (compress all of the output using ZLIB)

Random numbers:

u <int> (random number seed value for the simulated genome, default=11)

v <int> (random number seed value for generating somatic mutations on the simulated gneome, default=149)

w <int> (random number seed value for generating reads from the simulated tumor genomes, default=97)

Reads:

-c <int> (coverage of the genome, default = 0 [>0 distables option -n])

-n <int> (number of read pairs, default = 1000000)

-1 <int> (length of the first read, default = 100)

-2 <int> (length of the second read, default = 100)

-a <float> (disgard if the fraction of ambiguous bases higher than #FLOAT, default = 0.05)

-i <int> (average insert size, i.e. outer distance, between the two ends, default = 500)

Somatic:

-E <float> (sequencing base error rate from simulated tumor genome, default = 0.02)

-s <float> (somatic mutation rate, default = 0.01)

-t <float> (ratio of homozygous variants (the rest are heterozygous), default = 0.3333)

-p <float> (probabiltiy of observing tumor reads at any somatic site [bionomial distribution], default = 0.9)

-q <float> (fraction of somatic indels, default=0.15)

-x <float> (probability an somatic indel is extended, default = 0.3)

-h <float> (homozygous rate for somatic variants, default = 0.3333)

Germline:

-e <float> (sequencing base error rate from simulated normal genome, default = 0.02)

-d <int> (standard deviation of insert size, default = 50)

-r <float> (rate of mutations, including substitutions and indels, default = 0.001)

-R <float> (fraction of indels, default = 0.15)

-X <float> (probability an indel is extended, default = 0.3)

-H <float> (homozygous variant rate, default = 0.3333)

 

 使い方はgsimと同様。

SNVSniffer ssim reference.fa reads_base -e 0.01 -c 30

 

 

 

 

参考

SNVSnifferの解説ページに、SNVSniffer以外のツールについても、normalとtumorを比較しtumor specificな変異を検出するコマンドが記載されている。覚書としてここに載せておきます。

準備 

正常組織。

#解説ページではbwa memでなくbwa alnでマッピングしている。
bwa aln -t 4 reference.fa normal_1.fa.gz >tmp1.sai
bwa aln -t 4 reference.fa normal_2.fq.gz >tmp2.sai

#saiからpaired.samを作り、bamに変換。
bwa sampe -r '@RG id:gatech SM:illumina' -a 700 reference.fa tmp1.sai tmp2.sai normal_1.fq.gz normal_2.fq.gz | samtools view -bhS -@ 10 - >normal.bam

#sortとindex
samtools sort -@ 10 normal.bam > normal_sorted.bam
samtools index normal_sorted.bam

#pileupは時間がかかる
samtools mpileup -s -C 50 -f reference.fa normal_sorted.bam >normal.pileup

腫瘍サンプルのマッピングとpileupも同様の手順で進める。

 

tumorの変異コール。

 Varscan

java -Xmx4g -jar VarScan.jar somatic normal.mpileup tumor.mpileup out.vcf -output-vcf 1 --strand-filter --tumor-purity 0.9 --min-coverage-normal 3 --min-coverage-tumor 3

JointSNVMix ( paper

java -Xmx4g -jar VarScan.jar somatic normal.mpileup tumor.mpileup out.vcf -output-vcf 1 --strand-filter --tumor-purity 0.9 --min-coverage-normal 3 --min-coverage-tumor 3

SomaticSniper (HP

bam-somaticsniper -F vcf -f reference.fa tumor.sorted.bam normal.sorted.bam out.vcf

Mutect

java -Xmx16g -jar muTect.jar --analysis_type MuTect --reference_sequence reference.fa --input_file:normal normal.sorted.bam --input_file:tumor tumor.sorted.bam --vcf out.vcf

SNVSniffer

SNVSniffer somatic normal_header.sam tumor_header.sam normal.sorted.pileup tumor.sorted.pileup -o out.vcf

 

引用

SNVSniffer: an integrated caller for germline and somatic single-nucleotide and indel mutations

Liu Y, Loewer M, Aluru S, Schmidt B

BMC Syst Biol. 2016 Aug 1;10 Suppl 2:47