VEP

-

VEP(variant effect predictor)というソフトウェアがある。使い方だけ簡単にメモ

vep \
  --cache \
  --dir ../../../bin/.vep \
  --i input.vcf.gz \
  --o output.vep \
  --offline \
  --force_overwrite

としておくととりあえずPTV等は検出してくれる。

コメント

コメントを投稿

このブログの人気の投稿

Inverse-normal transformation

-
Inverse normal transformationとは正規分布しないデータをランク変換してそのランクを正規分布の累積分布関数で変換すること。変換後の値は正規分布する。Rでやると qnorm((rank(x)-0.5)/length(x)) となる。正規分布しない量的形質に正規性を持たせるためにこのような処理を行う。 BMIを目的とし、年齢・性別・バッチでで補正してQTLをしたい場合、まずlmで残差を計算して、この残差をInverse normal transformationして形質として回帰分析を行う。 res = lm(BMI〜Age+Sex+Batch)$residuals std_res = qnorm((rank(res)-0.5)/length(res)) のような感じ。 バラ色の解決策でもない様子 https://www.ncbi.nlm.nih.gov/pubmed/19526352
Previous »

SKAT

-
RのSKATパッケージを使いたい。 Z：行にバリアント、列にサンプルを持つ行列とする。値は0/1/2。 X：共変量行列 まず、null modelを作り(obj)これにテストしたいバリアント行列(Z)を投げる。目的変数が連続値の場合 out_type="C" を指定し、二値変数の場合 out_type="B" を指定する。 obj = SKAT_Null_Model(y.c ~ X, out_type="C") SKAT(Z, obj)$p.value obj = SKAT_Null_Model(y.b ~ X, out_type="C") SKAT(Z, obj)$p.value 上記のモデルはサンプル数が少ない時にやや保守的であるらしい。最近のバージョンではパーミュテーションp値を返すSKATBinaryという関数が実装されている。 out = SKATbinary(Z, obj) SKAT(Z, obj)$p.value SNPごとに効果量を変えるカーネルを使いたい場合、 kernel で指定する。 SKAT(Z, obj, kernel = "linear.weighted")$p.value デフォルトでは beta(maf, 1, 25) に従う重みを与えている。もし、 beta(maf, 0.5, 0.5) に従う重みを与えたい場合は以下のようにして変更できる。ちなみにこれはかなりゆっくりとした勾配になる。 SKAT(Z, obj, kernel = "linear.weighted", weights.beta = c(0.5, 0.5))$p.value SKAT-O SKAT-Oとは、SKATとburden-testをブレンドした検定方法である。ブレンド率 ρ \rho ρ は r.corr でコントロールする。 r.corr=0 のとき: SKATに等しい r.corr=1 のとき; Burden testに等しい method="SKATO" を実行すると ρ \rho ρ を色々に変えて検定を行い、一番小さなp値を返す(とマニュアルにあるが、実際実験してみるとそうならな...
Previous »

smartPCA

-
pupulation geneticsではPCAをよく使う。外れサンプルの検出や集団構造化の補正など。使い方メモ。 まず、Plink形式のファイルからEIGENSTRAT形式に変換しないといけない。それぞれのファイルの対応は以下の通り。 hoge.bed = fuga.geno hoge.fam = fuga.ind hoge.bim = fuga.snp これを変換するソフト convertf というソフトが同梱されていて、以下のような変換ファイルを食べさせて変換させることができる。 bed2eigen.par genotypename: hoge.bed snpname: hoge.bim indivname: hoge.fam outputformat: EIGENSTRAT genotypeoutname: fuga.geno snpoutname: fuga.snp indivoutname: fuga.ind familynames: NO convertf -p bed2eigen.par こうするとfuga.geno/fuga.snp/fuga.indができる。次にPCAするときは同様に.parファイルを作る。 eigen2pca.par genotypename: fuga.geno snpname: fuga.snp indivname: fuga.ind evecoutname: piyo.evec evaloutname: piyo.eval smartpca -p eigen2pca.par こうすると各サンプルの固有値の記載された piyo.eval と固有ベクトルの記載された piyo.evec が得られる。 さらにプロットするperlのスクリプトまでついており、 perl ploteig \ -i piyo.evec \ -c 1:2 \ -s ../out/myplot \ -x で ../out/myplot.pdf などが出力される。オプションは bed2eigen.par や eigen2pca.par ファイルに記載すれば良い。fuga.indファイルに人種情報を入れておくと自動で色分けしてくれたりして楽しいが通常のfamファイルは人種情報を保持できないため後から自分でindフ...
Previous »