joinに失敗

-
https://qiita.com/taruhachi/items/04ccca5e3bc20ad01791　より
joinが嫌いだ。並べども並べども join: hoge: is not sorted:とかいうから。これはLANGの設定の問題で、sort時とjoin時の言語設定により生じる誤差を行っているらしい。強制的にLANG=Cを指定することで回避できる。
LANG=C sort hoge > sorted_hoge
LANG=C sort fuga > sorted_fuga
LANG=C join sorted_hoge sorted_fuga > hogefuga
すみません勉強不足でした。

コメント

コメントを投稿

このブログの人気の投稿

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 »

locuszoom

-
使いたい。まず、マーカー名と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 すげー遅い
Previous »