kuroの覚え書き

96の個人的覚え書き

Nicotiana benthamiana

イネやアラビはゲノムデータが充実していて定法通りでNGSデータをマッピングできるのだけど、それ以外の植物だとどうだろう。
タバコ属植物Nicotiana benthamianaの場合、ゲノムがきちんと解析完了しているわけではなく、いまだcontig情報の塊でしかない。
transcriptのシーケンスデータとゲノムのショットガンシーケンスデータがあるのだが、RNAseqから遺伝子発現解析をするにはどういうストラテジーが有効か?
一つはtranscriptのfastaをreferenceとして、RNAseqのデータをマッピングしていく方法。ただし、この場合、transcriptデータが本当に全ゲノム情報をカバーしている保証はない。
もう一つはゲノムのcontigをreferenceにマッピングを実施し、gene modelのアノテーションファイルで遺伝子(transcript)を同定する方法。こちらの場合は発現領域をcontigがカバーしきれていない可能性もある。
簡単なのは前者なのだろうけど、情報量は後者のほうが多そうだ。

まずはreferenceを集めるところから。
https://solgenomics.net
ここから辿ろうとするとどうもリンクが切れてしまっていてたどり着けない。
https://btiscience.org/our-research/research-facilities/research-resources/nicotiana-benthamiana/
こちらからftp://ftp.solgenomics.net/genomes/Nicotiana_benthamiana/
ここに接続してファイルを貰ってくることにする。
transcriptのデータは
Annotation/Niben101/Niben101_annotation.transcripts.fasta.gz
genome contigは
Assemblies/Niben.genome.v1.0.1.contigs.fasta.gz
を利用。transcriptに対するgene modelのアノテーションはあるわけがなく、contigに対しては
Annotation/Niben101/Niben101_annotation.gene_models.gff
を使うことになりそうだ。

ところでcufflinksのオプションにアノテーションファイルを渡すとき、gtf形式のファイルを渡すように書籍などでは書かれている。
今回見つかったファイルはgffという拡張子になっていて、別物ではあるのだけれど、大体似た書式で、相互に変換も可能だということなので、その方法を探る。(どうもgffのままでも良さそうな感じもあるのだが)
これにはperlの変換スクリプトが広く使われているそうなので、それをもらってきた。
いろいろな人がいろいろに作成されているようなのだが、ひとまず次のものを試してみる。
http://seqanswers.com/forums/showpost.php?p=22529&postcount=4
アノテーションのソースだけちょいと書き換えて走らせようとしたら、perlのモジュールが無いよと言われたので、インストールすることに。
https://bi.biopapyrus.jp/perl/bioperl/
こちらを参考にBioPerlをインストール。

単純に並べ替えをするだけでなく、何やら抽出も行っているようで、出来上がったファイルサイズはかなり小さくなった。果たしてこれでいいのか?
というかgffとgtfでは最初の8カラムは全く同じで、9カラム目だけが書式が異なっているだけのようで、cufflinksでは9カラム目を使わない?


いっぽう、IGV.jsでBAMを表示するときにreferenceとして使うBEDファイルと言うものがあるが、コレも他の書式のファイルから生成することが可能なようなので、やってみる。
その1 fa.faiファイルから生成fastaから samtools faidxでまずindexを作成した上で

$ awk '{print $1 "\t0\t" $2}' input.fa.fai > output.bed

単にfaiから2カラム抜いてきてtab 0 tabで区切っているだけのようだ。いいのかそれで?

また、gtfからbedを作成するときは

$ cat input.gtf | awk '{OFS = "\t"} {print $1,$4,$5,$3,$6,$7}' > output.bed

faiから作るときは3カラムなのにgtf/gffから作るときは6カラム。コレはいかに?