bwaでtranscript.faにマッピングしたあとcufflinksで解析するには

bwaを使ってRNAseqのリードをtranscript.faをリファレンスとしてマッピング、出来上がったBAMファイルからcufflinksでFPKMを得ようとしたが、gtfファイルが無いので、FPKMを計算すると、マッピング領域のゆらぎを反映して1つのtranscriptの中にいくつかのlocusを勝手に想定してFPKMが計算されてしまうということがおこる。

これを防いでtranscriptごとにFPKMを計算させるにはやはりgtfファイルを用意するしかなさそうだ。
入り口はtranscripts.faファイル。
まずは

$ samtools faidx transcripts.fa

でfastaのindexを作成する。
transcriptsのfastaなのでヘッダ行にtranscript名、次の行に配列、となっているのでidxファイルは
transcript名　配列の長さ　配列の改行数
という形になっているので、これをgtfのフォーマットに合わせてawkで整形してやる。

$ awk '{print $1 "\tmarker\tCDS\t1\t" $2 "\t.\t+\t.\ttranscript_id \"" $1 "\"; gene_id \"" $1 "\";"}' transcripts.fa.fai > transcripts.gtf

gene_id, transcript_idをどちらもヘッダ行の内容にしておく。
出来上がったgtfファイルを用いてcufflinksを実行する。