投稿

8月, 2018の投稿を表示しています

shellでtsvをiterateする

タブ区切りでいくつかの要素のあるファイルを一行ずつ読み込み、カラムごとに別々の役割を持たせて処理をしたい。shellでの配列を使うとうまくいく。具体的には sample.txt という chr id dist pos 22 rs100000 0 1000000 21 rs200000 0 2000000 みたいなどこかで見たようなファイルがあった時、chrで参照先ファイルを選んで、その参照先ファイル内でidを検索したいというタスクがあるとする。参照先ファイルはchr22.txtみたいな形になっているとする。一行ずつアレイで読み込んで処理することになる。 while read line; do arr=(`echo line`) grep {arr[1]} chr{arr[0]}.txt done < sample.txt {arr[0]} もしくは arr とすると一番目の要素、 {arr[1]} とすると二番目の要素、となる。ちなみにアレイの内容すべてを表示するのは arr[*] もしくは arr[@] である。

vimで選択範囲をヤンクしたテキストで置き換える

vimで選択した範囲をテキストで置き換えたいんだ。 ctrl-v + y でテキストをヤンクして再度 ctrl-v で選択した範囲と置換したいが、したいのだが、選択してから思う。 x や d ではバッファが置き換えられてしまうと・・・。 みんな思うことのようであり ここ に"Stamping"として書いてあった。結論は vnoremap S "_d"0P である。 S はもともと”行削除して入力モード”だが Di でいいから置き換える。stampingと覚えよう。快適。

teeで中間ファイルをやっつける

ファイルを一行ずつ読み込んで 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

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

Rsq

某ソフトウェアで頻用される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オーバーフロー

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

vcfのINFOのparse

ワンライナー 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 たぶん遅い。