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

insectcell.exblog.jp
ブログトップ

タグ:シェル ( 5 ) タグの人気記事


2014年 05月 23日

genome.fa とtranscripts.gtf でtranscripts.fa をつくる


genome.fa とtranscripts.gtf でtranscripts.fa をつくる

編集 | 削除
tophat-cufflinks を使ってtranscripts.gtf は手に入れたもののtranscripts.fa がないから遺伝子見るの大変、とか、何でか知らないけどtranscripts.gtf しか貰えなかったパターンのときにtranscripts.fa を作らないといけない状況になります。

#transcripts.gtf とgenome.fa からtranscripts.fa を作る。
#transcripts.gtf を加工1
more transcripts.gtf | awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$9}' | grep exon | awk '{print $1"\t"$3"\t"$4"\t"$6"\t"$8"\t"$10}' | tr -d "\"" | tr -d ";" | awk '{print "transcript_id_"$5"\t"$1"\t"$2"\t"$3"\t"$6}' | awk '{print $1"\t"$2"\tsubstr($2,"$3","$4-$3+1")\t"$5/100"0"}' | awk '{print $1"\t"$2"\t"$3"\tn"substr($4,3,2)}' | sed -e '/n01/s/subs/@subs/g' | awk '{print $3}' | tr -d "\n" | tr "@" "\n" | grep substr > temp1.txt
#transcripts.gtf の加工2
more transcripts.gtf | awk -F "\t" '{print $1"\t"$3"\t"$4"\t"$5"\t"$9}' | grep exon | awk '{print $1"\t"$3"\t"$4"\t"$6"\t"$8"\t"$10}' | tr -d "\"" | tr -d ";" | awk '{print "transcript_id_"$5"\t"$1"\t"$2"\t"$3"\t"$6}' | awk '{print $1"\t"$2}' | uniq > temp2.txt

#genome.fa を1行にしておく
more genome.fa | sed -e '/>/s/$/@/g' | tr -d "\n" | tr ">" "\n" > line_genome.txt

#transcripts.gtf の情報を使って
paste temp2.txt temp1.txt | head -n 10
transcript_id_CUFF.1.1NODE_1000_length_89123_cov_46.018040substr($2,2,2497)
transcript_id_CUFF.2.1NODE_1000_length_89123_cov_46.018040substr($2,3291,1003)substr($2,4409,646)

#こういうコマンドをつくりたい
more line_genome.txt | grep NODE_1000_length_89123_cov_46.018040 | awk -F "@" '{print ">transcript_id_CUFF.1.1@"substr($2,2,2497)}' | tr "@" "\n" >> transcripts.fa
>transcript_id_CUFF.1.1
CGTAGTAGTATCGATCGTATCTACTATCGTACGTAGCTAGTCGATGCTAGTCGTAGCTAACTCGATCGATCGTAGTCTGTTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATCGTAGCTGATCGTAGCTGATGCTGTAGCTGATCGTAGCTGATCGTAGCTAGTCTAGCTGATGCTAGCTGATCGTAGCTAGCTGATCGTAGCTAGCTGATGCTAGCTGATCGTAGCTGATCGTAGCTAGCTGATCgTAGCTAGCTTATCGTGTAGTCGTAGCTAGCTGATCGTAGCTAGATAAAAACGGGGATCGATCGTAGCTAAAAAGCTTTTTTTGAAAAATCGTAGCTGATGCTAGCTGATGCTGTAGTGTGTGATTTATATTACTAGCTAGCCCCATGCTAGTCGATCACACAC

#transcripts.gtf の情報を使ってコマンドを作る
paste temp2.txt temp1.txt | awk '{print "more line_genome.txt | grep "$2" | awk -F \"@\" %{print \">"$1"@\""$3"}% | tr \"@\" \"\\n\" >> transcripts.fa"}' | tr "%" "'" > make_transcripts.sh

#準備ができたら実行
sh make_transcripts.sh

#できる
ls

genome.fa
transcripts.gtf
temp1.txt
temp2.txt
line_genome.txt
make_transcripts.sh
transcripts.fa <- NEW!!!


終わってから気づいたけどこれでいいじゃん

http://cufflinks.cbcb.umd.edu/gff.html

Extracting transcript sequences
The gffread utility can be used to generate a FASTA file with the DNA sequences for all transcripts in a GFF file. For this operation a fasta file with the genomic sequences have to be provided as well. For example, one might want to extract the sequence of all transfrags assembled from a Cufflinks assembly session. This can be accomplished with a command line like this:

gffread -w transcripts.fa -g /path/to/genome.fa transcripts.gtf

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

[PR]

by koretoki | 2014-05-23 15:00
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

[PR]

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


セル内改行は手と目で殺した
[PR]

by koretoki | 2013-12-18 13:31
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