前提：ゲノム配列の詳細な解析は実施されている。遺伝子予測はCDS領域に限定して実施されている。
やりたいこと：5’３’UTRがどこまでか推測し、cDNA領域のアノテーションを作る。
1. Trinityでde novoアセンブリを実行:
Trinityを使用してRNA-seqデータからde novoアセンブリを行う。
Trinity --CPU 8 --max_memory 24G --seqType fq --left reads_1.fq --right reads_2.fq --output trinity_output
single end readのときは
Trinity --CPU 8 --max_memory 24G --seqType fq --single reads_1.fq,reads_2.fq --output trinity_output
複数のfastqファイルを使うときはコンマで区切る必要がある。reads_*.fqのようにすると間がスペースとなるのでエラーが出るようだ。
2. TransDecoderを使用してCDS（Coding Sequence）を抽出:
Trinityの出力からCDSを抽出するために、TransDecoderを使用する。
TransDecoder.LongOrfs -t trinity_output/Trinity.fasta TransDecoder.Predict -t trinity_output/Trinity.fasta
これにより、CDSが Trinity.fasta.transdecoder.pep と Trinity.fasta.transdecoder.cds に保存される。
3. BLASTを使用して Trinity.fasta.transdecoder.cdsアセンブリを参照ゲノムにマッピング:
Trinityのアセンブリを参照ゲノムにマッピングし、遺伝子の対応付けを行う。
blastn -num_threads 8 -query Trinity.fasta.transdecoder.cds -db reference_genome -outfmt 6 -out trinity_vs_reference.blast
4. Blast結果のBED形式への変換:
blast2bedツールを使用してBlast結果をBED形式に変換し、BEDToolsを利用してUTR領域を抽出する。
blast2bed < trinity_vs_reference.blast > trinity_vs_reference.bed bedtools intersect -a trinity_vs_reference.bed -b reference_annotation.gff -wa -wb > trinity_vs_reference_annotated.bed
5. シェルスクリプトを使用してGFFファイルを生成:
シェルスクリプトを使用して、BEDファイルからGFFファイルを生成する。
以下は、この目的のための簡単なスクリプトの例。
awk 'BEGIN {OFS="\t"} {print $4, "custom_UTR_annotation", "UTR", $2+1, $3, ".", $6, ".", "gene_id " $1 "; transcript_id " $1}' trinity_vs_reference_annotated.bed > custom_UTR_annotation.gffこれでどうかな。
さらにCDSとUTRを一つにまとめたcDNAアノテーションにするには。
6. CDS領域のアノテーションGFF3とUTR領域のアノテーションGFF3の結合:
cat reference_annotation.gff custom_UTR_annotation.gff > merged_annotation.gff3
7. GFF3ファイルとして整形
awk -F'\t' 'BEGIN{OFS=FS} {gsub(/Parent=[^;]+;/,"",$9); if($3=="CDS"){print $1,$2,"exon",$4,$5,$6,$7,$8,$9} else print}' merged_annotation.gff3 > cdna_annotation.gff3まだ試行錯誤中。
このままではうまくいかない。
