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

insectcell.exblog.jp
ブログトップ
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


<< awkとsedで塩基配列をアミ...      内部寄生蜂のvenom をLC... >>