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

insectcell.exblog.jp
ブログトップ
2013年 06月 24日

awkとsedで塩基配列をアミノ酸配列にする

手元のマルチファスタな塩基配列をアミノ酸配列に変えたい。
でも、色んなパッケージを取ってくるのはよくわかんないし依存関係とか考えてるより培養したい...

awk とかsed とか誰のmac にでも入ってるものだけでなんとかしてみました。
どんな環境でも動きますけど死ぬほど遅いです。
ひとつの配列を処理するのに1秒ぐらいかかります。

使う時は、#########に囲われてる部分をコピーして"nucl2prot.sh" という名前で保存。
手元のマルチファスタな塩基配列を"transcripts.fasta" という名前にして

sh nucl2prot.sh transcripts.fasta proteins.fasta

と打てば、"proteins.fasta" という名前のマルチファスタなアミノ酸配列が手に入ります。

################################################################################################

#nucl2prot.sh

#塩基配列をアミノ酸に。コドンテーブルを目的に合わせて書き直せ
#数字の3桁"," 区切りを利用しているので合わない時は考えて
#すごく遅い
#$1=[input.file] $2=[output.file]

more $1 | grep -v ">" | tr "ATCGN" "12345" | sed -e :a -e 's/\(.*[0-9]\)\([0-9]\{3\}\)/\1,\2/;ta' | tr "12345" "ATCGN" | awk '{print $1","}' | sed -e 's/TCA,/S/g' -e 's/TCC,/S/g' -e 's/TCG,/S/g' -e 's/TCT,/S/g' -e 's/TTC,/F/g' -e 's/TTT,/F/g' -e 's/TTA,/L/g' -e 's/TTG,/L/g' -e 's/TAC,/Y/g' -e 's/TAT,/Y/g' -e 's/TAA,/*/g' -e 's/TAG,/*/g' -e 's/TGC,/C/g' -e 's/TGT,/C/g' -e 's/TGA,/*/g' -e 's/TGG,/W/g' -e 's/CTA,/L/g' -e 's/CTC,/L/g' -e 's/CTG,/L/g' -e 's/CTT,/L/g' -e 's/CCA,/P/g' -e 's/CCC,/P/g' -e 's/CCG,/P/g' -e 's/CCT,/P/g' -e 's/CAC,/H/g' -e 's/CAT,/H/g' -e 's/CAA,/Q/g' -e 's/CAG,/Q/g' -e 's/CGA,/R/g' -e 's/CGC,/R/g' -e 's/CGG,/R/g' -e 's/CGT,/R/g' -e 's/ATA,/I/g' -e 's/ATC,/I/g' -e 's/ATT,/I/g' -e 's/ATG,/M/g' -e 's/ACA,/T/g' -e 's/ACC,/T/g' -e 's/ACG,/T/g' -e 's/ACT,/T/g' -e 's/AAC,/N/g' -e 's/AAT,/N/g' -e 's/AAA,/K/g' -e 's/AAG,/K/g' -e 's/AGC,/S/g' -e 's/AGT,/S/g' -e 's/AGA,/R/g' -e 's/AGG,/R/g' -e 's/GTA,/V/g' -e 's/GTC,/V/g' -e 's/GTG,/V/g' -e 's/GTT,/V/g' -e 's/GCA,/A/g' -e 's/GCC,/A/g' -e 's/GCG,/A/g' -e 's/GCT,/A/g' -e 's/GAC,/D/g' -e 's/GAT,/D/g' -e 's/GAA,/E/g' -e 's/GAG,/E/g' -e 's/GGA,/G/g' -e 's/GGC,/G/g' -e 's/GGG,/G/g' -e 's/GGT,/G/g' -e 's/---,/-/g' -e 's/NNN,/N/g' -e 's/N..,/N/g' -e 's/.N.,/N/g' -e 's/..N,/N/g' > temp_nucl2prot_seq.txt

more $1 | grep ">" > temp_nucl2prot_head.txt

paste temp_nucl2prot_head.txt temp_nucl2prot_seq.txt | awk '{print $1"\n"$2}' > $2

rm temp_nucl2prot_head.txt
rm temp_nucl2prot_seq.txt

echo "(^q^)/"


################################################################################################

以上です。
ほんとに遅いから心配になったら

wc -l temp_nucl2prot_seq.txt

とかやって進行状況を見ても安心して下さい。

#20140109追記
塩基配列の文字数が3の倍数になってないとダメです。
その場合は生成されたアミノ酸配列に「,」が入ってきます。
数字の桁区切りの「,」を利用して3文字を区切っているため、例えば4文字であれば
ATCG→1234→1,234→A,S
と残念な結果になります。
5文字だと、
AATCG→11234→11,234→AA,S

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

[PR]

# by koretoki | 2013-06-24 21:31
2013年 06月 23日

ゲノムアノテーションファイルがエクセルだった

いただいたゲノムアノテーションファイルがエクセルでした。
genome.faとそれでtranscripts.faなファイルを作りたいと思います。
以下作業記録

# カラムとして管理していく。IDは付けない。
# まずエクセルから欲しそうなところを取り出してhenkan.txt にした。

# henkan.txt はこういうファイルでタブ区切り。エクセルから取ってきた。
# $1は遺伝子の名前、$2はアミノ酸配列、$3はどっちの鎖に付けられるか、$4ポジションのカラムでスタートとエンドが..で繋がってる。ポジションで複数のエキソンから成る場合は半角スペースで区切られてる。$5は全体のスタート位置、$6は全体のエンド位置を示す。$7はエキソンの数、$7は$4と同じ。$8はアノテーションが付いてる
Gene#Size (aa)StrandPositionsstartendexon#PositionsDescription
scaffold00001.g1.t1306-297..121429712141297..1214-
scaffold00001.g10.t1427+28319..295992831929599128319..29599-
scaffold00001.g100.t1364-224230..2253212242302253211224230..225321-
scaffold00001.g1000.t1337-2042669..20436792042669204367912042669..2043679-
scaffold00001.g1001.t11174+2043981..20475022043981204750212043981..2047502hypothetical protein hogehoge
scaffold00001.g1002.t11106+2047818..2048726 2048852..2049681 2049836..20514142047818205141432047818..2048726 2048852..2049681 2049836..2051414genehogehoge
scaffold00001.g173.t1914-371907..373954 374592..375274 375331..3753413719073753413371907..373954 374592..375274 375331..375341Putative protein of unknown function [hoge hoge]
scaffold00001.g18.t1645-43949..44175 44371..44416 44513..461744394946174343949..44175 44371..44416 44513..46174hogehoge

# 変換元のデータから場所を示すデータを切り出す。
more henkan.txt | grep scaffold | awk -F "\t" '{print $4}' | tr " " "." | awk -F "." '{print $1".."$3"."$4+0".."$6+0"."$7+0".."$9+0"."$10+0".."$12+0"."$13+0".."$15+0"."$16+0".."$18+0"."$19+0".."$21+0}' > kiridashi.txt


# kiridashi.txt の中身。
# 最大で7つのエキソンがあるので、startとendを7セット指定している(s1)..(e1).(s2)..(e2).(s3)..(e3).(s4)..(e4).(s5)..(e5).(s6)..(e6).(s7)..(e7)
head -n 2 kiridashi.txt
297..1214.0..0.0..0.0..0.0..0.0..0.0..0
28319..29599.0..0.0..0.0..0.0..0.0..0.0..0

# 何番染色体かを示す
head -n 2 scaffold_number.txt
1
1

# 上記の2ファイルを結合させて、stard end の指定と染色体番号を示すファイルを作る。
paste kiridashi.txt scaffold_number.txt | tr "\t" "." > temp002.txt

head -n 2 temp002.txt
297..1214.0..0.0..0.0..0.0..0.0..0.0..0.1
28319..29599.0..0.0..0.0..0.0..0.0..0.0..0.1

# ゲノムファイルを変換して一行のファイルにしてある。もちろん疑似配列な
more temp003.txt | awk '{print substr($1,1,50)}'
scaffold00001@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00002@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00003@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00004@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00005@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00006@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00007@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00008@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00010@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00011@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00012@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00013@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...
scaffold00014@CATTCTATTCTATAAACATTTTTTTTTTTTTCCGTA...


# 上のファイルを使って、こういうコマンドラインをつくりたい。
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,297,918)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,28319,1281)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt
more temp003.txt | grep 1 | awk -F "@" '{print substr($2,224230,1092)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)substr($2,0,0)}' >> transcriptSEQ.txt

# 上記のコマンドラインを作るためのコマンドライン
# 植木算的な意味で+1してる。エキソンが7未満の奴はawk のところだと($2,0,1)になってしまうので、sedでごまかした。
more temp002.txt | awk -F "." '{print "more temp003.txt | grep "$22" | awk -F \"@\" _{print substr($2,"$1","$3-$1+1")substr($2,"$4","$6-$4+1")substr($2,"$7","$9-$7+1")substr($2,"$10","$12-$10+1")substr($2,"$13","$15-$13+1")substr($2,"$16","$18-$16+1")substr($2,"$19","$21-$19+1")}_ >> transcriptSEQ.txt"}' | tr "_" "'" | sed -e 's/($2,0,1)/($2,0,0)/g' > make_transcriptsSEQ.sh

#-鎖用のファイルを作るよ
tr "ATCG" "TAGC" < transcriptSEQ.txt | rev > minus_transcriptSEQ.txt

#プラス鎖かマイナス鎖かを指定するファイル
more henkan.txt | grep scaffold | awk -F "\t" '{print $3}' > plus_or_minus.txt

#中身
more henkan.txt | grep scaffold | awk -F "\t" '{print $3}' | head -n 3
-
+
-

# プラスとマイナスを見てどっちを使うのか指定する。
paste plus_or_minus.txt transcriptSEQ.txt minus_transcriptSEQ.txt | awk '{if ($1=="+"){print $2} else if ($1=="-"){print $3}}' > transcriptSEQUENCE.txt

#遺伝子の名前取ってくる
more henkan.txt | grep scaffold | awk -F "\t" '{print $1}' | tr "." "_" > gene_names.txt

#合体させて欲しかったtranscripts.fasta をつくる
paste gene_names.txt transcriptSEQUENCE.txt | awk '{print ">"$1"\n"$2}' > transcripts.fasta


これで完成です。
ちゃんときれいにやる方法を勉強しないとなーと思い続けて2年目...



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

[PR]

# by koretoki | 2013-06-23 18:20
2013年 06月 15日

内部寄生蜂のvenom をLC-MSとRNA-seqで調べたよ

内部寄生蜂というのは、宿主の体内に住んで育つ蜂のグループのことです。
当然宿主の体内に入ってれば、宿主の免疫系から「おまえだれ?」って聞かれたりするわけで、侵入してるのがバレたらボコボコにされてしまいます。

内部寄生蜂の親は子どもが宿主の体内でボコられないようにvenom (なんか"毒"とは直訳しない雰囲気ある)を宿主の体内に打ち込んで、宿主の血球によう包囲を抑制しようとします。

キイロショウジョウバエを宿主として寄生する蜂がいて、キイロショウジョウバエはモデル生物で情報が多いので、キイロショウジョウバエ側の免疫の研究は結構されてました。でも、寄生蜂の方はほとんど情報がなかったので、今回著者らはvenom がどんなタンパク質なのかを調べました。

PLoS One. 2013 May 23;8(5):e64125. doi: 10.1371/journal.pone.0064125. Print 2013.
Integrative approach reveals composition of endoparasitoid wasp venoms.
Goecks J, Mortimer NT, Mobley JA, Bowersock GJ, Taylor J, Schlenke TA.


実験では2種類の寄生蜂、L. b. とL. h. を使いました。
RNA-seq で寄生蜂腹部で発現してる遺伝子を特定し、LC-MS で特定したvenomu 線の中身のペプチドをそこにはりつけることで、venom のタンパク質を特定しようとしています。

RNA-seq から作ったトランスクリプトデータに、mass で得たアミノ酸配列を貼付けたところ、L.b. で129のvenom 遺伝子候補、L.h. で176のvenom 遺伝子候補がそれぞれとれました。

それらの中から、雌蜂特異的に発現するもの(雄vs雌のRT-PCR比較から)、
venom 腺特異的に発現するもの(雌体vsvenom 腺のRT-PCR比較から)を絞り込むことにしました。
※他の蜂のvenomとして、Calreticulin や、Heat-shock protein 70 が知られているので、上記の方法ではそれらの遺伝子のような、他の組織で他の仕事をしているvenom 遺伝子を取りこぼすリスクがあるそうです。

ランダムに10個のvenom 遺伝子候補を選んでRT-PCRしたところ、
L.b. でもL.h. でも10個中の8個は雌特異的に発現をしていました。

L.b. ではそれら8個の遺伝子はvenom 腺以外のところでも発現していました。L.b. では、ショウジョウバエ寄生蜂Asobara japonica と同様に、venom 遺伝子が卵巣など他の組織でも発現しているのでしょう。
一方で、L.h. では8個の遺伝子はvenom 腺特異的な発現をしていました。


次に、既知のvenom 遺伝子と相同性のある遺伝子を集めて、分泌シグナルついてるやつだけを集めました。
GO解析をしたところ、L.b. のサンプルでは、antioxidant activity がエンリッチしていました。
ショウジョウバエの免疫応答では、酸素ラジカルを発生させていると考えられているので、L.b. は抗酸化タンパクをvenom として用いて、免疫に対抗していると考えられます。



RNA-seq とMASS を組み合わせた方法は、ゲノムが分かっていない生物の分泌タンパクを特定するのにすごい効くと著者らは述べています。




寄生蜂の話をするといっつもvenom っていう謎物質が出てきて苦しめられてきたのですが、蜂によって全然違うんですね。
今回の研究で見つかった200個くらいのvenom 候補をひとつづつ試してもらって、venom の成分完全同定されるようになると期待しています。
あと、いままでGO階席でぱりっと意味をとる論文に出会えてなかったので、すごい絞り込んだ遺伝子群でG0解析→エンリッチみっけ→エンリッチさせている遺伝子これな→知見と会わせてばっちり考察してるのが楽しかったです。





以下メモ

L.b., L. h. について
幼虫か蛹に産卵
L. b. は卵を免疫から隠して、血球が卵の表面にふれるのを防ぐ。宿主の免疫は寄生されていても強い。
L. h. は宿主の免疫を破壊する。

最初はcDNAをクローニングして進めた。4つとれた。でも、100くらいはありそうなので、ほんの表面だけしか見れていない。

P. puparum のmass もやったけどトランスクリプトがないから56個の予測venom のうち12個しか取れなかった。

Chelonus inanitus のESTsやった。29個みつかって、いくつかは完全にユニークな配列だった。

Microctonus hyperodae でRNA-seqやった。8個の全長venom 遺伝子見つかった。

Nasonia vitripennis でmass やってゲノムのにはっつけたら79venom 遺伝子見つかった。

→核酸+アミノ酸のアプローチがいい
→NGSでせめるのがいい



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

[PR]

# by koretoki | 2013-06-15 10:30
2013年 06月 13日

ショウジョウバエの脂肪酸合成酵素群の配列をまるっと欲しい場合

特定の機能に関連した遺伝子の配列がまるっと欲しい気持ちになることがしばしばあります。
そういうときはKEGGを使うのがいいよって教えていただいたのでときどき使っています。
例えば、ショウジョウバエの脂肪酸合成酵素群をまるっと欲しいときは、

KEGG pathway maps のメニュー(以下リンク)を開く
http://www.genome.jp/kegg-bin/get_htext?query=00061&htext=br08901.keg&option=-a

脂肪酸合成パスウェイが見たいので"fatty acid biosynthesis"をクリックして開く(以下リンク)
http://www.genome.jp/kegg-bin/show_pathway?map00061

"Reference patheway" のところをクリックしてショウジョウバエの選んでG0をクリック(以下リンク)
http://www.genome.jp/kegg-bin/show_pathway?org_name=dme&mapno=00061&mapscale=&show_description=hide

脂肪酸合成パスウェイのうち、ショウジョウバエで見つかっている物は緑色に塗られています。
緑に塗られているところをクリックするとそれらの配列を見ることができるので、それをポチポチコピペしていけば脂肪酸合成パスウェイ上の配列を集めることができます。



でもやりたくないです。
30個くらいやるのに1時間かかったことあるし。
twitter で聞いたら@nonshu さんが教えて下さったのでまとめました。


マップの番号が00061であることは分かっています。(マップの左下に書いてある)
キイロショウジョウバエの略称はdmeなので(マップ上のどれか適当に選んでクリックして配列のファスタのヘッダーに書いてあるのがそれ。例えば、">dme:Dmel_CG11198 ACC; Acetyl-CoA carboxylase (EC:6.3.4.14 6.4.1.2); K11262 acetyl-CoA carboxylase / biotin carboxylase [EC:6.4.1.2 6.3.4.14] (A) "。)

基本URL1[http://www.genome.jp/dbget-bin/get_linkdb?-t+genes+path:]のおしりに[dme00061]をくっつけてURLを書きます。

例)http://www.genome.jp/dbget-bin/get_linkdb?-t+genes+path:dme00061

すると、IDとDefinition がでてくるのでIDだけコピペして+で繋げます。
例)dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

それを、基本URL2[http://www.genome.jp/dbget-bin/www_bget?-f+-n+a+]のおしりにくっつける。

例)http://www.genome.jp/dbget-bin/www_bget?-f+-n+a+dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

リンク先に欲しい配列がファスタフォーマットで待ってます。

もし塩基で欲しい時には、aをnに変えます。
基本URL3[http://www.genome.jp/dbget-bin/www_bget?-f+-n+n+]のおしりにIDをくっつける。

例)http://www.genome.jp/dbget-bin/www_bget?-f+-n+n+dme:Dmel_CG11198+dme:Dmel_CG12170+dme:Dmel_CG17374+dme:Dmel_CG3523+dme:Dmel_CG7842

リンク先に欲しい配列がファスタ形式で待ってます。


はやく迎えに行ってあげて





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

[PR]

# by koretoki | 2013-06-13 19:52
2013年 02月 14日

経費な本リスト

1 2013/12/06 数量化革命 アルフレッドWクロスビー amazon
2 2013/12/12 ベイズ統計データ解析 姜興起 amazon
[PR]

# by koretoki | 2013-02-14 17:46