2014年 01月 09日
makeblastdb -in ref.fa -dbtype nucl -hash_index blastn -db ref.fa -query que.fa -outfmt 6 > blastout.txt makeblastdb -in ref.pep.fa -dbtype prot blastp -db ref.pep.fa -query que.pep.fa -outfmt 6 > blastout.txt #
by koretoki
| 2014-01-09 12:25
2014年 01月 07日
メタボローム解析の反復データをひとつにまとめたい。 n=3のものとn=1のものを比べたいのだけど、できればn=3全部使いたい。 (a1+a2+a3)/3 ってやっていいのかわからなかったから調べてみた。 #もともと持っているのはこんなデータ。いわゆるMartixとして得られたデータである。 name a1 a2 a3 CYT001 0.0007330321 0.0004886212 0.0008590029 CYT002 0.1763202000 0.1228157000 0.1730612000 CYT003 0.0137760700 0.0123575800 0.0164345100 CYT004 0.0121694300 0.0103101300 0.0095795930 CYT005 0.0008732506 0.0006902750 0.0008904565 CYT006 0.0007769990 0.0005277728 0.0005936377 CYT007 0.0010050660 0.0005964599 0.0010336850 CYT008 0.0000000000 0.0000000000 0.0000000000 CYT009 0.0000000000 0.0000000000 0.0000000000 CYT010 0.0005402510 0.0004378932 0.0003912269 CYT011 0.0000000000 0.0000000000 0.0000000000 CYT012 0.0000000000 0.0000000000 0.0000000000 #(これがずっと続いてる) #R様に読んでいただく。 > data<-read.table("data.txt",header=T,row.names=1) > attach(data) > summary(data) a1 a2 a3 Min. :0.0000000 Min. :0.0000000 Min. :0.0000000 1st Qu.:0.0003296 1st Qu.:0.0002691 1st Qu.:0.0000955 Median :0.0015751 Median :0.0015310 Median :0.0010649 Mean :0.0696770 Mean :0.0660012 Mean :0.0445601 3rd Qu.:0.0200691 3rd Qu.:0.0198329 3rd Qu.:0.0128484 Max. :2.4163290 Max. :2.2341410 Max. :1.5474730 #読んでいただけた。実際にはデータが全部表示されるが省略。 > data a1 a2 a3 CYT001 0.0007330321 0.0004886212 0.0008590029 CYT002 0.1763202000 0.1228157000 0.1730612000 CYT003 0.0137760700 0.0123575800 0.0164345100 CYT004 0.0121694300 0.0103101300 0.0095795930 CYT005 0.0008732506 0.0006902750 0.0008904565 CYT006 0.0007769990 0.0005277728 0.0005936377 CYT007 0.0010050660 0.0005964599 0.0010336850 #総和がこんなにズレてる... > sum(a1) [1] 15.53798 > sum(a2) [1] 14.71827 > sum(a3) [1] 9.936896 #これをそのままプロットして1:1のところに線引いてみる。 > plot(a1,a3) > abline(0,1,col="red") #ずれずれっぽいけどよく見えないからlog > plot(log(a1),log(a3)) > abline(0,1,col="red") #とりあえず総和で割ればいけそう > plot(log(a1/sum(a1)),log(a3/sum(a3))) > abline(0,1,col="red") (a1/sum(a1)+a2/sum(a2)+a3/sum(a3))/3 すればとりあえず物質間の比については使えるデータになりそう。 a1*x+a2*y+a3*zの形にしたい。 総和でも良さそうだけど怖いから回帰で比を出してもらうことにした。 #それぞれ回帰してもらう #回帰モデルする > fit_12<-lm(a1~a2) > fit_23<-lm(a2~a3) #中見る1vs2 > fit_12 Call: lm(formula = a1 ~ a2) Coefficients: (Intercept) a2 -0.0005735 1.0643829 #中見る2vs3 > fit_23 Call: lm(formula = a2 ~ a3) Coefficients: (Intercept) a3 0.006123 1.343758 #調節する時にかける値は x:y=1:1.0643829 y:z=1:1.343758 x:y:z=1:1.0643829:1.430273=0.2861512:0.3045745:0.4092743 0.2861512*a1+0.3045745*a2+0.4092743*a3 で欲しい値が取れる。 #総和の比でやると > sum(a1) [1] 15.53798 > sum(a2) [1] 14.71827 > sum(a3) [1] 9.936896 > sum(a1)/sum(a2) [1] 1.055693 > sum(a2)/sum(a3) [1] 1.481174 > #x:y=1:1.055693 > #y:z=1:1.481174 > 1.055693*1.481174 [1] 1.563665 > #x:y:z=1:1.055693:1.563665 > 1*1/(1+1.055693+1.563665) [1] 0.2762921 > 1*1.055693/(1+1.055693+1.563665) [1] 0.2916796 > 1*1.563665/(1+1.055693+1.563665) [1] 0.4320283 > #x:y:z=1:1.055693:1.563665=0.2762921:0.2916796:0.4320283 総和の方が強い?補正になってるのは、値の大きな子が引っ張ってるのかなー? #回帰で求めた比がちゃんとハマるか心配 > plot(a1,a2) > abline(0,1/1.0643829,col="red") > plot(a2,a3) > abline(0,1/1.343758,col="red") > lm(a3~a1) Call: lm(formula = a3 ~ a1) Coefficients: (Intercept) a1 -0.00182 0.66564 > plot(a3,a1) > abline(0,1/0.66564,col="red") #logでもハマるか一応 > plot(log(a1),log(a2)) > e<-exp(1) > e [1] 2.718282 > abline(-log(1.0643829),1,col="red") > plot(log(a2),log(a3)) > abline(-log(1.343758),1,col="red") > abline(0,1,col="blue") > plot(log(a3),log(a1)) > abline(0,1,col="blue") > abline(-log(0.66564),1,col="red") 赤線が回帰の線 青線は1:1の線 人気ブログランキングに参加しています。 応援よろしくお願いします。 #
by koretoki
| 2014-01-07 01:18
2013年 12月 20日
手元にあるマルチファスタファイルをシングルファスタのたくさんのファイルにします。 今回使うマルチファスタファイルは multi.fasta という名前のファイルです。 それをバラバラにして、遺伝子名.fasta という名前のファイルにします。 #使うものつくる $ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" > m2s.sh #動くようにする $ chmod 777 m2s.sh #動かす $ sh m2s.sh 使ったもの:more, grep, tr, awk, sed, chmod, sh ##################################################### #multi.fasta の中を見てみる $ more multi.fasta >gene_a CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC >gene_b TGACTGATGCTATACGACGATGCTAGCTACGTAGCTAGCTAGCTAGCTGATGAATATATTTATTTCTGCTAGCTATAGCTGATGCTAGTCGTAGTGGGGGCTAGCTAGTCGTAGCTAGTCGTAGTCGTAGCTGATGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA >gene_c CGATGCTATTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCAAAAAAAAGCTGTAGCTAGCATCTATCTACTATCGTATCACTGACGGCGGCGCGGCGCGCGGGCCGGAGAGCGGACGTCATACGTCAGTACTGTCAGTCAGTACTGCAGTACGTCGTGTCGTAGTCATGACTGTACGTGTAGCACGTCGTACGTACGTCGAGCGAGCTGATTCGATCTGATCTGATTCGTAGCGGCGTGCGCGATCGTAGCTATCGATGCTGATCGAGCTGAACTGGTACTCTGGTAGCTAGCTAGCTTGACGTCATGTACGTCATGTCAGTCAGTACTGTCAGTACGTACTGACTGACGGTCAGTACG #gene_a のみを切り出してシングルファスタファイルにしたかったらこう $ more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_a | awk -F "@" '{print ">"$1"\n"$2}' > gene_a.fasta #gene_a.fasta をみてみる $ more gene_a.fasta >gene_a CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC #これを全部の遺伝子について書こうとするとこう $ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" | more more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_a | awk -F "@" '{print ">"$1"\n"$2}' > gene_a.fasta more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_b | awk -F "@" '{print ">"$1"\n"$2}' > gene_b.fasta more multi.fasta | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" | grep gene_c | awk -F "@" '{print ">"$1"\n"$2}' > gene_c.fasta #それをm2s.sh という名前で保存します $ more multi.fasta | grep ">" | tr -d ">" | awk '{print "more multi.fasta | sed -e %/>/s/$/@/g% | tr -d \"\\n\" | tr \"\>\" \"\\n\" | grep "$1" | awk -F \"@\" %{print \"\>\"$1\"\\n\"$2}% > "$1".fasta"}' | tr "%" "'" > m2s.sh #m2s.sh を実行可能ファイルにします $ chmod 777 m2s.sh #ファイルが揃いました $ ls gene_a.fastam2s.shmulti.fasta #さっきつくったgene_a.fasta いらない $ rm gene_a.fasta #すっきり $ ls m2s.shmulti.fasta #m2s.sh の実行 $ sh m2s.sh #実行結果。遺伝子毎にバラバラになってる $ ls gene_a.fastagene_b.fastagene_c.fastam2s.shmulti.fasta #中身を確認 $ head gene_* ==> gene_a.fasta <== >gene_a CTAGCTAGTCTAGCTATCTATCGTATCGTAGCTGATGCTAGTCGTAGCTGATCGTAGCTGATGCTGTAGCGATGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTGATCGTAGCTAGTCGATGTGTGTGTGGGGGCGCGCCGGCGCGGCGC ==> gene_b.fasta <== >gene_b TGACTGATGCTATACGACGATGCTAGCTACGTAGCTAGCTAGCTAGCTGATGAATATATTTATTTCTGCTAGCTATAGCTGATGCTAGTCGTAGTGGGGGCTAGCTAGTCGTAGCTAGTCGTAGTCGTAGCTGATGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA ==> gene_c.fasta <== >gene_c CGATGCTATTTTTTTTTTTTTTTTTTTCCCCCCCCCCCCCAAAAAAAAGCTGTAGCTAGCATCTATCTACTATCGTATCACTGACGGCGGCGCGGCGCGCGGGCCGGAGAGCGGACGTCATACGTCAGTACTGTCAGTCAGTACTGCAGTACGTCGTGTCGTAGTCATGACTGTACGTGTAGCACGTCGTACGTACGTCGAGCGAGCTGATTCGATCTGATCTGATTCGTAGCGGCGTGCGCGATCGTAGCTATCGATGCTGATCGAGCTGAACTGGTACTCTGGTAGCTAGCTAGCTTGACGTCATGTACGTCATGTCAGTCAGTACTGTCAGTACGTACTGACTGACGGTCAGTACG 実際のファスタファイルにはスペースが入っていたりすることが多いので、事前に殺すか-F "\t" で回避。 元のマルチファスタと同じところにシングルファスタがもりもり入るとうざいから事前に mkdir single として箱つくって > "$1".fasta" のところを > single/"$1".fasta" にするとすっきり 人気ブログランキングに参加しています。 応援よろしくお願いします。 #
by koretoki
| 2013-12-20 14:24
2013年 12月 18日
#同じ値に複数の化合物が含まれているデータのモデルデータ $ more test2 C00001 ID_1 ID_2 ID_3 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04 C011010 ID_1 ID_2 ID_3 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D. C001983,C03231,C04213 ID_1 ID_2 ID_3 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02 C034143 ID_1 ID_2 ID_3 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D. C044555,C022233 ID_1 ID_2 ID_3 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D. #ひとまずいらないIDを消して化合物毎にわけます。 $ more test2 | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' | tr "," "\n" | awk '{print NR"\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7}' | sort -nr 8 C022233 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D. 7 C044555 6 C034143 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D. 5 C04213 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02 4 C03231 3 C001983 2 C011010 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D. 1 C00001 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04 #データが空になっちゃった化合物のところにデータを戻します。 $ more test2 | awk '{print $1"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10}' | tr "," "\n" | awk '{print NR"\t"$1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6"\t"$7}' | sort -nr | awk '(NR==1){a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8; print; next}{if($3<0){print $1"\t"$2"\t"a} if($3>0){print; a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8} }' | more 8 C022233 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D. 7 C044555 4.8 2.1E-04 4.8E-04 1.3E-04 N.D. N.D. 6 C034143 6.8E-10 8.1E-99 N.D. N.D. N.D. N.D. 5 C04213 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02 4 C03231 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02 3 C001983 N.D. 2.1E-05 1.1E-02 0.9E-02 1.1E-02 1.1E-02 2 C011010 1.11 6.4E-04 6.3E-04 N.D. N.D. N.D. 1 C00001 2.1E-04 10.1E-04 0.4E-04 N.D. N.D. 1.2E-04 #前の行を参照する方法 awk ' (NR==1){a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8; print; next} #1行目は絶対入ってるのでaに値を代入してそのまま吐いて次の行にいく { if($3<0){print $1"\t"$2"\t"a} #行のデータが空だったらaを吐く if($3>0){print; a=$3"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8} #行のデータが空じゃなかったらそのまま吐いて自分の行の値をa に代入する }' | more セル内改行は手と目で殺した #
by koretoki
| 2013-12-18 13:31
2013年 12月 12日
genome.fasta->geneprediction.gff->codingseq.fasta export AUGUSTUS_CONFIG_PATH=~/tools/augustus/augustus-2.4/config/ ~/tools/augustus/augustus-2.4/src/augustus --cds=on --codingseq=on --species=homo_n0rr homo_n0rr.fasta > homo_n0rr.gff more homo_n0rr.gff | grep "#" | tr -d "#" | tr -d " " | sed -e 's/startgene/startgene@/g' -e 's/endgene/endgene@/g' | tr -d "\n" | tr "@" "\n" | grep endgene | tr "]" "[" | sed -e 's/coding/[coding/g' | tr atcg ATCG | awk -F "[" '{print ">augustus_hn_"$1"\n"$3}' > codingseq_homo_n0rr.fasta #$5でpepseq 人気ブログランキングに参加しています。 応援よろしくお願いします。 #
by koretoki
| 2013-12-12 01:13
|
アバウト
以前の記事
2021年 02月 2018年 02月 2017年 12月 2016年 12月 2016年 06月 2014年 09月 2014年 05月 2014年 01月 2013年 12月 2013年 06月 2013年 02月 2012年 11月 2011年 07月 2011年 05月 2010年 11月 2010年 05月 2010年 04月 2009年 12月 2009年 09月 2009年 08月 2009年 07月 2009年 06月 2009年 04月 2009年 03月 2009年 02月 2009年 01月 2008年 12月 2008年 11月 カテゴリ
最新のトラックバック
タグ
その他のジャンル
ブログパーツ
ファン
記事ランキング
ブログジャンル
画像一覧
|
ファン申請 |
||