位相/スイッチエラー率を計算するための検索ツール

Wouter De Coster

2019-03-13 02:53:51 UTC

view on stackexchange narkive permalink

真理vcfファイルとテストvcfファイルを指定して、位相/スイッチエラー率を計算するツールを探しています。 WhatsHapを使用してvcfのフェーズを実行し、その結果を、私が持っているグラウンドトゥルースフェーズのvcfと比較したいと思います。ツールが見つからないので、自分で書きたくありません。かなり一般的なことのようです...それは存在すると思いますが、今日のgoogle-fuは弱いです。

テスト行の例：

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samplechr10100588。 GA。 PASS AF = 0.18; AC = 1; NS = 150; AN = 2; EAS_AF = 0.14; EUR_AF = 0.39; AFR_AF = 0.08; AMR_AF = 0.31; SAS_AF = 0.2; VT = SNP; DP = 11907 GT：PS 1 | 0 ：48005chr10102385。 AG。 PASS AF = 0.07; AC = 1; NS = 150; AN = 2; EAS_AF = 0.02; EUR_AF = 0.0; AFR_AF = 0.1; AMR_AF = 0.19; SAS_AF = 0.01; VT = SNP; DP = 22665 GT：PS 0 | 1 ：48005chr10105170。 GA。 PASS AF = 0.03; AC = 1; NS = 150; AN = 2; EAS_AF = 0.02; EUR_AF = 0.0; AFR_AF = 0.0; AMR_AF = 0.19; SAS_AF = 0.0; VT = SNP; DP = 17140 GT：PS 0 | 1 ：48005chr10105365。 GT。 PASS AF = 0.19; AC = 1; NS = 150; AN = 2; EAS_AF = 0.14; EUR_AF = 0.39; AFR_AF = 0.09; AMR_AF = 0.31; SAS_AF = 0.2; VT = SNP; DP = 21845 GT：PS 1 | 0 ：48005chr10106057。 TC。 PASS AF = 0.13; AC = 1; NS = 150; AN = 2; EAS_AF = 0.02; EUR_AF = 0.0; AFR_AF = 0.32; AMR_AF = 0.19; SAS_AF = 0.01; VT = SNP; DP = 22559 GT：PS 0 | 1 ：48005

真実の行の例

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samplechr10100554。 C CT、CTT30。 。 GT：AD 1 | 2：0,1,1chr10100588。 G A30。 。 GT：AD 0 | 1：1,1chr10102385。 A G30。 。 GT：AD 1 | 0：1,1chr10102636。 T C30。 。 GT：AD 1 | 1：0,2chr10102757。 T C30。 。 GT：AD 0 | 1：1,1
chr10104815。 GCC G、GC30。 。 GT：AD 2 | 1：0,1,1chr10105170。 G A30。 。 GT：AD 1 | 0：1,1chr10105365。 G T30。 。 GT：AD 0 | 1：1,1chr10106057。 T C30。 。 GT：AD 1 | 0：1,1chr10106110。 C G30。 。 GT：AD 1 | 1：0,2chr10106261。 A G30。 。 GT：AD 1 | 0：1,1chr10108612。 T G30。 。 GT：AD 1 | 0：1,1chr10108646。 AT30。 。 GT：AD 0 | 1：1,1chr10108834。 G A30。 。 GT：AD 1 | 0：1,1chr10110840。 C T30。 。 GT：AD 1 | 1：0,2chr10111743。 AT30。 。 GT：AD 1 | 0：1,1chr10112016。 T A30。 。 GT：AD 1 | 0：1,1chr10112262。 C T30。 。 GT：AD 1 | 1：0,2chr10113006。 C T30。 。 GT：AD 0 | 1：1,1chr10113031。 G A30。 。 GT：AD 0 | 1：1,1chr10113136。 G A30。 。 GT：AD 1 | 0：1,1chr10113359。 A T30。 。 GT：AD 1 | 1：0,2chr10113583。 C G30。 。 GT：AD 0 | 1：1,1chr10113997。 G C30。 。 GT：AD 1 | 1：0,2

テストする各ファイルの数行を教えてください。 2つのvcfファイルの同じバリアントの遺伝子型フィールドを比較したいだけですか？

質問を編集して、例の行を含めました。トゥルースファイルには、すべてのバリアントが段階的に含まれています（合成二倍体メソッドhttps://www.nature.com/articles/s41592-018-0054-7およびhttps://github.com/lh3/CHM-evalに基づいて作成されています） / tree / master / dip-call）。

$ awk -F "\ t" -vOFS = "\ t" 'BEGIN {print "#CHROM"、 "POS"、 "ID"、 "REF"、 " ALT "、" Truth "、" Test "} {if（/ ^ [^＃] /）{var = $ 1 $ 2 $ 4 $ 5; gt = substr（$ 10,1,3）; if（NR == FNR）{a [var] = gt; next;} if（var in && a [var]！= gt）{print $ 1、$ 2、$ 3、$ 4、$ span> 5、gt、a [var]}}} 'truth.vcf test.vcf #CHROM POS ID REF ALT Truth Testchr10100588。 G A 1 | 0 0 | 1chr10102385。 A G 0 | 1 1 | 0chr10105170。 G A 0 | 1 1 | 0chr10105365。 G T 1 | 0 0 | 1chr10106057。 TC 0 | 1 1 | 0

awk -F "\ t" -vOFS = "\ t" 'BEGIN {print "#CHROM"、 "POS"、 "ID"、 "REF"、 "ALT"、 "Truth"、 "Test"} {if（/ ^ [^＃] /）{var = $ 1 $ span> 2 $ 4 $ span> 5; gt = substr（ $ 10,1,3）; if（NR == FNR）{a [var] = gt;次; } if（var in && a [var]！= gt）{print $ 1、$ 2、$ 3、$ 4、$ span> 5、gt、a [var]}}} 'truth.vcf test.vcf

awk -F "\ t" -vOFS = "\ t" 'BEGIN {print "#CHROM"、 "POS"、 "ID" 、 "REF"、 "ALT"、 "Truth"、 "Test"} {if（/ ^ [^＃] /）{var = $ 1 $ span> 2 $ 4 $ span> 5; gt = substr（ $ 10,1,3）; if（NR == FNR）{a [var] = gt;次; } tot ++; if（var in && a [var]！= gt）{bad ++; print $ 1、$ 2、$ 3、$ 4、$ span> 5、gt、a [var]}}} END {printf "正解：％s of％s（％s％）\ n"、tot-bad、tot 、（tot-bad）* 100 / tot; printf "Wrong：％s of％s（％s％）\ n"、bad、tot、bad * 100 / tot;} 'truth.vcf test.vcf

#CHROM POS ID REF ALT Truth Testchr10100588。 G A 1 | 0 0 | 1chr10102385。 A G 0 | 1 1 | 0chr10105170。 G A 0 | 1 1 | 0chr10105365。 G T 1 | 0 0 | 1chr10106057。 T C 0 | 1 1 | 0正解：0/5（0％）不正解：5/5（100％）