人気ブログランキング | 話題のタグを見る

したっぱ昆虫細胞研究者のメモ

insectcell.exblog.jp
ブログトップ
2014年 01月 09日

makeblastdb

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")
メタボデータをくっつけた_e0160319_1135217.png

#ずれずれっぽいけどよく見えないからlog
> plot(log(a1),log(a3))
> abline(0,1,col="red")
メタボデータをくっつけた_e0160319_114586.png


#とりあえず総和で割ればいけそう
> plot(log(a1/sum(a1)),log(a3/sum(a3)))
> abline(0,1,col="red")
メタボデータをくっつけた_e0160319_1152369.png


(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")

メタボデータをくっつけた_e0160319_116455.png


赤線が回帰の線
青線は1:1の線

人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking


# 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" にするとすっきり



人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking


# by koretoki | 2013-12-20 14:24
2013年 12月 18日

エクセルなメタボロームデータを化合物IDで揃えたテーブルにしたい

#同じ値に複数の化合物が含まれているデータのモデルデータ
$ 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日

gen. pre.

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

人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking


# by koretoki | 2013-12-12 01:13