mach2qtl

minimacでimputationしたdosageのgenotypeファイルでassociationしたい。とくに量的変数。

このためにmach2qtlというソフトをminimacの開発者がリリースしているが、使い方に癖があるのでメモ

–datfile sample.dat

C age
C female
T lg10CRP

Cで共変量、Tで目的変数を指定

–pedfile sample.ped

# FID	IID	Father Mother	sex	age	female	lg10CRP
4001	4001	0	0	m	41.591	0	0.580
4014	4014	0	0	m	61.864	0	1.021
4021	4021	0	0	m	53.859	0	0.845

PLINKでいうところのfamファイルである。6列目以下に目的変数もしくは共変量を記載する

–dosefile sample.dose.vcf.gz

これが現行のminimac3の出力ファイルを受けつけてくれない。DosageConvererをつかって形式を変更する。

dosageConvertor \
  --vcfDose sample.dose.vcf.gz \
  --info sample.info \
  --prefix sample_mldose \
  --type mach \
  --format 1

–dosefile sample_mldose.dose.vcf.gz

としてようやく実行

mach2qtl \
  --datfile sample.dat \
  --pedfile sample.ped \
  --infofile sample.info \
  --dosefile sample_mldose.dose.vcf.gz \
  --rsqcutoff 0.7 > output.txt

つまり標準出力で返る。リダイレクトで受けるのだが、フォーマットがやや利用しにくい。きれいなテーブルにしたいなら下記のスクリプトを使用してください。

cat sample_mldose.doseassoc | \
awk '{
  if(NR==1){print "TRAIT\tCHR\tPOS\tALLELES\tFREQ1\tRSQR\tEFFECT2\tSTDERR\tCHISQ\tPVALUE\tN"}
  if($1~/^lg10CRP$/ && $2~/:/) {
    split($2, arr, ":");
    print($1"\t"arr[1]"\t"arr[2]"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10)
  }
}' > sample_mldose.doseassoc.txt

lg10CRPは自分のフェノタイプ名に変更すること。