$ head Niben.genome.v1.0.1.contigs.fasta >Niben101Scf00001Ctg001 AAAAAAAGGATTAAGTGTCATAAATGTGGTAAATTTGGTCATTATGCAAGTGAGTGTAAA ACTCAGGAAAATATTAAGAGTCTAGATTTAGATGATAAACTTAAGGATTCTTTGTGTAAG ATTCTACTAAATTCTGATTATAGTTCTGATGTATCTGATTCCTCTTCTACTGGTTCTATT TCTGATGAAAATCTAAGGGTTTTACATGATGAAGATTTATATTCTGATTCTGACAGATCT TTTTCGAATGATGATTGTGCTTGCAAGAAAGANGATGAAGATTTATATTCTGATTCTGAC AGATCTTTTTCGAATGATGATTGTGCTTGCAAGAAAGATGATAGTTTGTATGATTTAATT TCTAATTTTCAGGATCTGAATGTGAACATGATTAATGATGAGAAGATAATAAAGATTTTA AAAACAGTAGAAGATAAAGATTTGAGAAATAAAATCTTGGATATTGTTTTGACTAATGAT GAGCCTAATACCAATTTTCCCTCTTCAAATTTTAGTTTAGAACAGCCTGATTACTATTCG
コンティグごとの配列が並んでいる。
gffファイルは
$ head Niben101_annotation.gene_models.gff Niben101Ctg00001 . contig 1 500 . . . ID=Niben101Ctg00001 Niben101Ctg00002 . contig 1 500 . . . ID=Niben101Ctg00002 Niben101Ctg00003 . contig 1 500 . . . ID=Niben101Ctg00003 Niben101Ctg00004 . contig 1 500 . . . ID=Niben101Ctg00004 Niben101Ctg00005 . contig 1 500 . . . ID=Niben101Ctg00005 Niben101Ctg00006 . contig 1 500 . . . ID=Niben101Ctg00006 Niben101Ctg00007 . contig 1 500 . . . ID=Niben101Ctg00007 Niben101Ctg00008 . contig 1 500 . . . ID=Niben101Ctg00008 Niben101Ctg00009 . contig 1 500 . . . ID=Niben101Ctg00009 Niben101Ctg00010 . contig 1 500 . . . ID=Niben101Ctg00010
gffから変換したgtfではcontigの情報とかは消されており
$ head Niben101_annotation.gene_models.gtf Niben101Ctg00054 . exon 167 487 . + . gene_id "Niben101Ctg00054g00001";transcript_id "Niben101Ctg00054g00001.1"; Niben101Ctg00074 . exon 60 113 . - . gene_id "Niben101Ctg00074g00004";transcript_id "Niben101Ctg00074g00004.1"; Niben101Ctg00074 . exon 186 495 . - . gene_id "Niben101Ctg00074g00004";transcript_id "Niben101Ctg00074g00004.1"; Niben101Ctg00116 . exon 2 25 . - . gene_id "Niben101Ctg00116g00002";transcript_id "Niben101Ctg00116g00002.1"; Niben101Ctg00116 . exon 120 314 . - . gene_id "Niben101Ctg00116g00002";transcript_id "Niben101Ctg00116g00002.1"; Niben101Ctg00116 . exon 391 501 . - . gene_id "Niben101Ctg00116g00002";transcript_id "Niben101Ctg00116g00002.1"; Niben101Ctg00129 . exon 25 327 . - . gene_id "Niben101Ctg00129g00001";transcript_id "Niben101Ctg00129g00001.1"; Niben101Ctg00141 . exon 60 273 . + . gene_id "Niben101Ctg00141g00002";transcript_id "Niben101Ctg00141g00002.1"; Niben101Ctg00141 . exon 320 495 . + . gene_id "Niben101Ctg00141g00002";transcript_id "Niben101Ctg00141g00002.1"; Niben101Ctg00174 . exon 62 487 . + . gene_id "Niben101Ctg00174g00001";transcript_id "Niben101Ctg00174g00001.1";
このようにexon情報だけになっている。
なおfaiを
$ samtools faidx Niben.genome.v1.0.1.contigs.fasta
と 生成して、中を確認すると
$ head Niben.genome.v1.0.1.contigs.fasta.fai Niben101Scf00001Ctg001 591 24 60 61 Niben101Scf00002Ctg001 301 649 60 61 Niben101Scf00002Ctg002 159 980 60 61 Niben101Scf00003Ctg001 744 1166 60 61 Niben101Scf00004Ctg001 207 1947 60 61 Niben101Scf00004Ctg002 750 2182 60 61 Niben101Scf00005Ctg001 43104 2969 60 61 Niben101Scf00005Ctg002 11715 46816 60 61 Niben101Scf00005Ctg003 8402 58751 60 61 Niben101Scf00005Ctg004 5331 67318 60 61 $ wc -l Niben.genome.v1.0.1.contigs.fasta.fai 288736 Niben.genome.v1.0.1.contigs.fasta.fai
このように288736のcontigが含まれている事がわかる。
一見、これでマッピングすればいけそうな気がするのだが、やってみるとどういうわけかアノテーションの付いたものはFPKMが全て0になってしまっている。なんでかなー。
アプリは突貫で一応できており、検索、BAMの表示もできるんだけどなあ
ちなみにgtfで抽出された行をgffで探してみると
Niben101Ctg00054 maker gene 167 487 . + . ID=Niben101Ctg00054g00001;Alias=snap_masked-Niben101Ctg00054-processed-gene-0.0 Niben101Ctg00054 maker mRNA 167 487 . + . ID=Niben101Ctg00054g00001.1;Parent=Niben101Ctg00054g00001;Alias=snap_masked-Niben101Ctg00054-processed-gene-0.0-mRNA-1;Note="Ethylene-responsive transcription factor 7";Ontology_term=GO:0003677,GO:0006355,GO:0003700 Niben101Ctg00054 maker exon 167 487 . + . ID=Niben101Ctg00054g00001.1:exon:001;Parent=Niben101Ctg00054g00001.1;Alias=snap_masked-Niben101Ctg00054-processed-gene-0.0-mRNA-1:exon:4812 Niben101Ctg00054 maker CDS 167 487 . + 0 ID=Niben101Ctg00054g00001.1:cds:001;Parent=Niben101Ctg00054g00001.1;Alias=snap_masked-Niben101Ctg00054-processed-gene-0.0-mRNA-1:cds
こんな感じになっている。
Niben101Ctg00054の167-487についてはgffではgene mRNA exon CDSの４つが乗っているが、変換したgtfではexonしか残っていないな。これがマッピングでアノテーションが全然反映されていない原因かな。
locusで検索をかけられるように細工してみるとたしかにマッピングはされているが、アノテーション情報は反映されていないようだ。というかよく見るとlocusの表記が一致していない。これが原因だね。
$ head Niben.genome.v1.0.1.scaffolds.nrcontigs.fasta.fai Niben101Scf00001 591 27 60 61 Niben101Scf00002 584 655 60 61 Niben101Scf00003 744 1276 60 61 Niben101Scf00004 1153 2060 60 61 Niben101Scf00005 274325 3260 60 61 Niben101Scf00006 349369 282185 60 61 Niben101Scf00007 604 637404 60 61 Niben101Scf00008 963 638046 60 61 Niben101Scf00009 790 639053 60 61 Niben101Scf00010 622004 639884 60 61
こっちだね。
あと、gff→gtfの変換だが
$ gffread Niben101_annotation.gene_models.gff -T -o Niben101_annotation.gene_models.gtf
とするとcufflinksに付属するツールで処理してくれる。perlとか要らない。
head Niben101_annotation.gene_models.gtf Niben101Ctg00054 maker exon 167 487 . + . transcript_id "Niben101Ctg00054g00001.1"; gene_id "Niben101Ctg00054g00001"; Niben101Ctg00054 maker CDS 167 487 . + 0 transcript_id "Niben101Ctg00054g00001.1"; gene_id "Niben101Ctg00054g00001"; Niben101Ctg00074 maker exon 60 113 . - . transcript_id "Niben101Ctg00074g00004.1"; gene_id "Niben101Ctg00074g00004"; Niben101Ctg00074 maker exon 186 495 . - . transcript_id "Niben101Ctg00074g00004.1"; gene_id "Niben101Ctg00074g00004"; Niben101Ctg00074 maker CDS 60 113 . - 0 transcript_id "Niben101Ctg00074g00004.1"; gene_id "Niben101Ctg00074g00004"; Niben101Ctg00074 maker CDS 186 287 . - 0 transcript_id "Niben101Ctg00074g00004.1"; gene_id "Niben101Ctg00074g00004"; Niben101Ctg00116 maker exon 2 25 . - . transcript_id "Niben101Ctg00116g00002.1"; gene_id "Niben101Ctg00116g00002"; Niben101Ctg00116 maker exon 120 314 . - . transcript_id "Niben101Ctg00116g00002.1"; gene_id "Niben101Ctg00116g00002"; Niben101Ctg00116 maker exon 391 501 . - . transcript_id "Niben101Ctg00116g00002.1"; gene_id "Niben101Ctg00116g00002"; Niben101Ctg00116 maker CDS 2 25 . - 0 transcript_id "Niben101Ctg00116g00002.1"; gene_id "Niben101Ctg00116g00002";
このようにcontigだけの行は除かれるが、CDSの情報は残るようだ。なにせcufflinksの純正ツールなんだから間違いないだろう。
さて、状況はすべて把握できた。いざマッピング。
バッチリじゃなかろうか。とりあえず出来上がったBAMをIGVでちゃんと表示できた。
なお、デスクトップ版のIGVを使うにあたって、Genomes＞Create .genome File...
で、上の
Niben.genome.v1.0.1.scaffolds.nrcontigs.fasta
Niben101_annotation.gene_models.gtf
をそれぞれFASTA file、Gene fileに指定してgenomeファイルを作成し、リファレンスに指定してやること。
