hanajori

プロセス置換というものを知った。無断な中間ファイルを作らなくてすむ。具体的には以下 並び直してjoinしたいとき join <(sort a.txt) <(sort b.txt) 何行あるかわからないコメントをのぞいてpasteしたいとき paste <(grep -v -E "^#" a.txt) <(grep -v -E "^#" b.txt) そのほかの活用例 例1

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値を返す(とマニュアルにあるが、実際実験してみるとそうならな...

ダイナミックリンクライブラリエラー再び。gcc-8.1.0を落としてきてlibstdc++をビルドし直す。以下、非ルートユーザでの実行例。 wget http://mirrors.concertpass.com/gcc/releases/gcc-8.1.0/gcc-8.1.0.tar.gz tar -xf gcc-8.1.0.tar.gz cd gcc-8.1.0 ./contrib/download_prerequisites mkdir build && cd build ../configure \ --enable-languages=c,c++ \ --prerix $HOME/local \ --disable-bootstrap \ --disable-multilib make -j 80 make install #環境変数で今ビルドしたライブラリを先にロードするように指示する export LD_PRELOAD=$HOME/local/lib64/libstdc++so.6

従来のShapeitやBeagleと違うアルゴリズムによる高速なhaplotype phasing software。サンプルが多いときに 有利 。 リファレンスなしでvcfをphasingする。リファレンスパネル作成のため。 eagle \ --vcfTarget=target.vcf.gz \ --numThreads=15 \ --chrom=9 \ --bpStart=100000 \ --bpEnd=200000 \ --geneticMapFile genetic_map_hg19.txt.gz \ --outPrefix=target.phased \ 2>&1 | tee example_ref.log リファレンスを用いてvcfをphasingする。imputation受け手のphasingなど。 eagle \ --vcfRef=reference.vcf.gz \ --vcfTarget=target.vcf.gz \ --numThreads=15 \ --chrom=9 \ --bpStart=100000 \ --bpEnd=200000 \ --geneticMapFile genetic_map_hg19.txt.gz \ --outPrefix=target.phased \ 2>&1 | tee example_ref.log

使いたい。まず、マーカー名とp-valueからなるmetal形式のインプットファイルを作る echo -e “MarkerName\tP-value” > …/tmp/GWAStutorial.metal tail -n +2 ../dat/assoc/GWAS.assoc | \ sed -E "s/ +/\t/g" | sed -E "s/^\t//" | \ awk '{print "chr"$1":"$3"\t"$9}' >> ../tmp/GWAS.metal input.metal Markername P-value rs11111 0.01 rs22222 0.02 chr1:123456 0.5 みたいにする。列名は必要で、これあれば指定なしに読み込んでくれる。rs noがなくても、 chr1:721290 の形式にしておけばin placeで表示してくれる。 base positionで領域を指定したい場合 locuszoom \ --metal input.metal --build hg19 --pop ASN --source 1000G_March2012 \ --chr chr9 --start 210000000 --end 230000000 Markernameで指定したい場合 locuszoom \ --metal input.metal --build hg19 --pop ASN --source 1000G_March2012 \ --refsnp rs12565286 --flank 500kb Genenameで指定したい場合 locuszoom \ --metal input.metal --build hg19 --pop ASN --source 1000G_March2012 \ --refgene ABCA1 --flank 500kb すげー遅い