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年目... 人気ブログランキングに参加しています。 応援よろしくお願いします。
by koretoki
| 2013-06-23 18:20
|
アバウト
以前の記事
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月 カテゴリ
最新のトラックバック
タグ
その他のジャンル
ブログパーツ
ファン
記事ランキング
ブログジャンル
画像一覧
|
ファン申請 |
||