hanajori

公式サイトのトップで遊べる 黄色い```>_```ボタンをクリックすると実行環境に入れる。matplotlibもインストールされていた（いらんけど）

手元端末(Mac)からsshで計算サーバ(RedHat)に繋いでいるが、そろそろパスワードを打つのが億劫になってきた。公開鍵認証をする。サーバのrootもあるし、 sshd_configをいじってだな、、、とか気合い入れたが、デフォルトで公開鍵は有効になってました。むしろRHEL6のマニュアルの該当箇所は「パスワード認証を禁止する」ところから始まっているという。 とりあえず手元で鍵を作る。 $ ssh-keygen -t rsa Enter file in which to save the key (/Users/user/.ssh/id_rsa): Enter passphrase (empty for no passphrase): Enter same passphrase again: と聞かれるので名前を変えたければ変える。変えなければ id_rsa と id_rsa.pub ができる。scpでもsftpでも何でもいいので id_rsa.pub をサーバに持って行って、 サーバの home/.ssh/authorized_keys に id_rsa.pub の内容を追記する。で サーバで パーミッションを変更する。 cat id_rsa.pub >> ~/.ssh/authorized_keys chmod 700 ~/.ssh chmod 600 ~/.ssh/authorized_keys これで ssh -i ~/.ssh/id_rsa user@setsuzoku-saki でパスワードなしで接続でき、2回目以降は ssh user@setsuzoku-saki で接続できるはずなのだが残念ながらうまくいかない。 ~/.ssh/config に Host setsuzoku HostName setsuzoku-saki User user IdentityFile ~/.ssh/id_rsa とするとよい。この場合は ssh setsuzoku みたいなので接続できる

流石に研究ログを取らねばならない。しかし、どうしたものか。ベストプラクティス（笑）で書いたような形でディレクトリを細切れに構成していくと、再現性は高いが可視性がよくない。かと言ってそれぞれのディレクトリの場所をテキストファイルにコピペしていくのはあまりにも生産性が低いように思う。せっかく考えてディレクトリを配置しているので。 それぞれのディレクトリに書いたREADME.mdを自動収集するスクリプトを試してみる。 getLogs () { PRJ=$1 PRJDIR=$2 # ディレクトリなければ作る mkdir -p /home/hanajori/logfiles/README_LOG/$PRJ/weekly mkdir -p /home/hanajori/logfiles/README_LOG/$PRJ/update 　#指定ディレクトリ以下のmdファイルを全て収集して表示、作成日時と共にweeekly.mdファイルに保存 for md in `find $PRJDIR -name README.md`; do ls -l$md | cut -f 6- -d " "; cat $md; printf "\n" >>$PRJDIR/README_LOG/weekly/`date +"${PJR}_%y%m%d"_weekly.md done 　cd$PRJDIR/README_LOG/weekly 　lastweek=`ls -lt $PRJDIR/README_LOG/weekly *weekly.md | head -2 | tail -1` 　thisweek=`ls -lt$PRJDIR/README_LOG/weekly *weekly.md | head -1 | tail -1` 　diff $lastweek$thisweek > ../update/`date + "${PRJ}_%y%m%d_update.md } #ログを撮りたいプロジェクトを記載 getLogs PROJECT_X /home/hanajori/PRJ_X getLogs PROJECT_Y /home...

JOBという名前のプロセスをkillしたい場合。一度topの出力をテキストに保存してからループで投げる top -b -d 1 -n 1 > ids.txt for id in `cat ids.txt | grep JOB | cut -d " " -f 1`; do kill $id done

vcfとくにvcf.gzから直接色々な解析をしてくれるソフトとして rvtest がある。 通常のassociation test(Wald test)を行いたい場合 rvtest \ --pheno example.pheno \ --inVcf example.vcf \ --single wald \ --out example.wald vcfは圧縮されていても良い。example.phenoはいわゆるplink.famファイル。6列目以降にもphenotypeデータをおくことができて、 --mpheno X で指定する。6列目がデフォルト、7列目が --mpheno 2 となる。column名で指定する場合には --pheno-name PHENOTYPE などとする。 rvtest \ --pheno example.pheno \ --inVcf example.vcf \ --single wald \ --covar example.covar \ --covar-name C1,C2 \ --mpheno 2 \ --out example.wald covariateを取り入れたい場合は --covar example.covar とする。 --covar-name C1,C2 などとしてcovariate名を指定する。 rvtest \ --pheno example.pheno \ --inVcf example.vcf \ --single wald \ --out example.wald Burden test(CMC)を実行する場合 ※tabixでindex化されたvcf.gzファイルが必要 rvtest \ --pheno example.pheno \ --invcf example.vcf.gz \ --setFile example.set \ --burden cmc \ --out example example.set ファイルを用意する。これは annotation1 1:1-3 annotation2 2:4-6 といった形式のファイルである。実行すると example.CMC.assoc というファイルが出力...

LSFというものを使っている。忘れないようにメモ書き 基本的に bsub -q QUE -o OUT "shell command" で動く。例えば bsub -q QUE -o OUT "pwd; ls" とすると、 OUT ファイルに pwd; ls の結果、つまりカレントディレクトリとlsの結果が出力される。処理したい内容を myscript.sh に書いておいて bsub -q mini_que -o my.log "bash myscript.sh" とすれば良い。たくさん投げる場合 for l in {1..22}; do bsub -q mini_que -o myout_$l.log "bash myscript.sh$l" done とかすれば良い。 gnu paralllel が使えれば bsub -q mini_que -o myout.log "parallel myscript.sh ::: {1..22}" が いけるよう だが、例によって入ってなくて使えない。

バイオインフォの解析などを共有したいと思っている。 バイオインフォ（にかぎらずデータ解析一般）では基本的に ・インプットデータ ・処理ソース ・アウトプットデータ が基本コンポーネントである。さらにいうとインプットから直接アウトプットにいけない場合も多いので、中間ファイルが必要な場面も多い。取り回しやデータの保全を考えるとインプットデータとアウトプットデータはそれぞれのディレクトリに単独で居てくれた方が嬉しいので、中間ファイルは別におかれるべきである。 考えた結果 projectA_1 |--dat # input data |--src # analysis source |--tmp # temporary file |--out # output data という配置でやっている。 そして、さらに下流の解析に行くときは projectA_2 |--dat |--src |--tmp |--out |--init.sh を作る。 projectA_1/out は projectA_2/dat の中身になる（ことが多い）。このときシムリンクにするのが良いか、コピーした方が良いかは考え中。コピーするメリットはデータが保全できること。リンクのメリットはディスクサイズとパイプラインとしての再利用性。このときのコピー（あるいはリンク貼り作業）も init.sh で記録しておくと後々再現性が高い。 src の中も処理内容によってスクリプトが別個になってしまうのは仕方がないとして、最終的に全ての処理を一括で行う main.sh みたいのがあるとよい。ファイル名に気をつけるのとエラー処理をきちんとするスクリプトを書いておけば、各ディレクトリの init.sh と main.sh を実行していけばパイプラインになる。あるいは誰かに渡すとき最悪これだけやってもらうだけで解析が進む。 この制約を自分に課すようになってから、自分のいま扱っているファイルが dat なのか tmp なのか out なのか意識するようになった（ような気がする）。制約は自由度を下げるが、逆に型にはめることでイレギュラーに敏感になることができる。トレーニングの間は有効に作用するこ...

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フ...