hanajori

ファイルを一行ずつ読み込んで sagasu_AがあればA.txt sagasu_BがあればB.txt へ流し込んでいきたい場合 cat oreno.txt | tee \ >(grep -E "sagasu_A" >> A.txt) \ >(grep -E "sagasu_B" >> B.txt) > /dev/null ちなみに、遅い

rmarkdownをrenderするのがめんどかった。さらに出来上がりを確認するのにわざわざhtmlを開くのもめんどかった。。。以下のスクリプトをパスの通ったディレクトリに置いておいて、 render oreno.rmd とすればrenderしてさらに勝手にブラウザでopenしてくれる。 !/bin/bash html=${1%%.rmd}.html Rscript -e "rmarkdown::render('$1')" open $html

某ソフトウェアで頻用されるRsqという値について v a r ( θ ) ρ ( 1 − ρ ) \frac{var(\theta)}{\rho(1-\rho)} ρ ( 1 − ρ ) v a r ( θ ) ​ と説明がある。ここで θ \theta θ はimputed dosageのベクトル、 ρ \rho ρ はminor allele frequencyである。この計算をRでやると Rsq = var(theta)/(rho*(1-rho)) とおもったらRのvarは不偏分散を使うので値が合わない（こら）。 Rsq = sum((theta - mean(theta))^2)/(rho*(1-rho)) ですね。

awkの数値比較は1e-308あたりでおかしくなるらしい。 https://lists.gnu.org/archive/html/bug-gawk/2015-04/msg00011.html 上記を参考にすると、1e-308を下回りそうな値については0を足すことで0と判定させることができ、この値と比較することで正しい数値比較ができる。 たとえば、ある値が1e-6以下なら表示させたいとする。しかし、この中に1e-308より小さいのがあると漏らす。 間違う例(1e-6以下なら表示したい) >> echo 1e-666 | awk '{if($1<1e-6) {print $1}}' 正しく判定できる例 >> echo 1e-666 | awk '{if(!($1+0)<1e-6) {print $1}}' 1e-666

ワンライナー extract () { cat $1 | awk -v regex="$2=([-|.|e|0-9]+);" \ '{match($8, regex, arr);\ if(length(arr)==0){print "NA"} else {print arr[1]}} ' } とかにするとよろしい。 for var in AC AF AN BaseQRankSum ; do extract oreno.vcf $var > $var.txt done たぶん遅い。

プロセス置換というものを知った。無断な中間ファイルを作らなくてすむ。具体的には以下 並び直して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値を返す(とマニュアルにあるが、実際実験してみるとそうならな...