「ほっ」と。キャンペーン
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


<< ダーウィンは一回だけevolv...      makeblastdb >>