Hatena::ブログ(Diary)

kuroの覚え書き

[Vmax] [Macintosh/PC] [Network] [Make] [science] [Programing]

注:個人的覚え書きであり、いかなる内容も保証されるものではありません。
(ツッコミは歓迎しますが、全てに応えられるとも限りませんのであしからず。)

2018-02-21

[][]データベースアプリ完成

いやあ、一つテンプレートとなるアプリを作っておけば、内容が違っていてもそれなりのデータベースアプリがすぐに作成できてしまうな。

前のやつには結局5ヶ月ほどかけたが、今回のはほぼ2週間で出来上がった。細部少々荒いし、マルチユーザーを意識せず、ローカル動作を前提にしていたりするので、ややこしい処理が減っているけどね。

さあ、解析環境の構築もできてきたので、実際にデータをあーでもないこーでもないといじるお楽しみの時間ですよ。

2018-02-16

[][] cummeRbundによるクラスター解析

以前に一通り試したクラスタリングを実際にやってみた。

忘れていたことも含めて覚書

> library("cummeRbund")    #ライブラリ読み込み
> cuff <- readCufflinks("cuffdiff_resultsのディレクトリ")    #cuffdiffの出力ディレクトリからデータを読み込み)
> setwd("~/")    #ワーキングディレクトリを設定
> my.genes <- Genes(cuff)
> mySigGeneIDs <- getSig(cuff,alpha=0.05,level="genes")    #p value<0.05のものを抽出。数が多い時はもっと絞ってみる。
> mySigGenes <- getGenes(cuff, mySigGeneIDs)
> ic <- csCluster(mySigGenes, k=8)    #8つのクラスタに分類
> icp <- csClusterPlot(ic)   #plotを作成
> icp    #plotを描画ー>pdf等で保存しておく
> old.op <- options(max.print=999999)    #textの出力リミットを外す
> sink("~/ic.txt")    #出力先を指定
> ic
> sink()

mySigGenesが多すぎるとクラスタリングに時間が掛かるし、結果もぼやける。結局上位の候補だけを解析することになるので、最初の抽出で数を絞るべきだ。

どうしても数が多いとき、icの結果を表示するとgetOption("max.print")を超えたぶんは表示されないし、データとして出力もできないので、old.op <- options(max.print=999999)でリミットをいっぱいまで上げてから出力する。

出力結果から該当するクラスタのメンバーをコピーしてexcel等で処理する。

2018-02-14

[][] cuffmergeがエラー

Error (GFaSeqGet): end coordinate (515) cannot be larger than sequence length 504
Error (GFaSeqGet): end coordinate (515) cannot be larger than sequence length 504
Error (GFaSeqGet): subsequence cannot be larger than 506
Error getting subseq for CUFF.13.1 (2..512)!
	[FAILED]
Error: could not execute cuffcompare

どういうわけかこんなエラーが出て失敗する。ちゃんとしたゲノムデータではなくてツギハギのコンティグデータをリファレンスにしているから、うまく動かないんだろうか?

[][]StringTieを入れる

TopHat→cufflinks→cuffdiff→cummeRbundというのがRNAseq解析の王道であるが、最近はHISAT2→StringTie→Ballgownというのが流行りはじめているらしい。これはやっとかないと。ということでまずは環境構築から。

StringTieのインストールはbrewでできるらしいので、早速

$ brew install stringtie

とやってみると、なにやらbrew自体がrubyで動いているのだけれど、rubyのバージョンが古いぞ、と怒られた。

/usr/local/Homebrew/Library/Homebrew/brew.rb:12:in `<main>': Homebrew must be run under Ruby 2.3! You're running 2.0.0. (RuntimeError)

で、OSXのrubyのアップデートをするならrbenvをbrewで入れて・・・っておい、そのbrewが動かないんだってば。どうすりゃいいのさ?

$ brew update

とりあえずコレでいいみたい。

せっかくなので

$ brew install rbenv
$ rbenv install 2.4.1
$ rbenv global 2.4.1

さて、気を取り直してstringtieを入れよう。

$ brew install stringtie
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 1 tap (homebrew/core).
==> New Formulae
fruit

Error: No available formula with the name "stringtie" 
==> Searching for a previously deleted formula (in the last month)...
Warning: homebrew/core is shallow clone. To get complete history run:
  git -C "$(brew --repo homebrew/core)" fetch --unshallow

Error: No previously deleted formula found.
==> Searching for similarly named formulae...
==> Searching local taps...
Error: No similarly named formulae found.
==> Searching taps...
==> Searching taps on GitHub...
Error: No formulae found in taps.

なんと。移転したらしい。移転先は教えてくれないのか・・・

brewsci/bioかbrewsci/scienceかな。

$ brew tap brewsci/science
==> Tapping brewsci/science
Cloning into '/usr/local/Homebrew/Library/Taps/brewsci/homebrew-science'...
remote: Counting objects: 598, done.
remote: Compressing objects: 100% (595/595), done.
remote: Total 598 (delta 1), reused 186 (delta 1), pack-reused 0
Receiving objects: 100% (598/598), 536.53 KiB | 681.00 KiB/s, done.
Resolving deltas: 100% (1/1), done.
Tapped 579 formulae (617 files, 1.6MB)
$ brew install stringtie
Updating Homebrew...
==> Installing stringtie from brewsci/science
==> Downloading https://homebrew.bintray.com/bottles-science/stringtie-1.3.3b.el_capitan.b
######################################################################## 100.0%
==> Pouring stringtie-1.3.3b.el_capitan.bottle.tar.gz
&#127866;  /usr/local/Cellar/stringtie/1.3.3b: 5 files, 509.8KB

こっちがアタリだった。

$ stringtie -h
StringTie v1.3.3b usage:
 stringtie <input.bam ..> [-G <guide_gff>] [-l <label>] [-o <out_gtf>] [-p <cpus>]
  [-v] [-a <min_anchor_len>] [-m <min_tlen>] [-j <min_anchor_cov>] [-f <min_iso>]
  [-C <coverage_file_name>] [-c <min_bundle_cov>] [-g <bdist>] [-u]
  [-e] [-x <seqid,..>] [-A <gene_abund.out>] [-h] {-B | -b <dir_path>} 
Assemble RNA-Seq alignments into potential transcripts.
 Options:
 --version : print just the version at stdout and exit
 -G reference annotation to use for guiding the assembly process (GTF/GFF3)
 --rf assume stranded library fr-firststrand
 --fr assume stranded library fr-secondstrand
 -l name prefix for output transcripts (default: STRG)
 -f minimum isoform fraction (default: 0.1)
 -m minimum assembled transcript length (default: 200)
 -o output path/file name for the assembled transcripts GTF (default: stdout)
 -a minimum anchor length for junctions (default: 10)
 -j minimum junction coverage (default: 1)
 -t disable trimming of predicted transcripts based on coverage
    (default: coverage trimming is enabled)
 -c minimum reads per bp coverage to consider for transcript assembly
    (default: 2.5)
 -v verbose (log bundle processing details)
 -g gap between read mappings triggering a new bundle (default: 50)
 -C output a file with reference transcripts that are covered by reads
 -M fraction of bundle allowed to be covered by multi-hit reads (default:0.95)
 -p number of threads (CPUs) to use (default: 1)
 -A gene abundance estimation output file
 -B enable output of Ballgown table files which will be created in the
    same directory as the output GTF (requires -G, -o recommended)
 -b enable output of Ballgown table files but these files will be 
    created under the directory path given as <dir_path>
 -e only estimate the abundance of given reference transcripts (requires -G)
 -x do not assemble any transcripts on the given reference sequence(s)
 -u no multi-mapping correction (default: correction enabled)
 -h print this usage message and exit

Transcript merge usage mode: 
  stringtie --merge [Options] { gtf_list | strg1.gtf ...}
With this option StringTie will assemble transcripts from multiple
input files generating a unified non-redundant set of isoforms. In this mode
the following options are available:
  -G <guide_gff>   reference annotation to include in the merging (GTF/GFF3)
  -o <out_gtf>     output file name for the merged transcripts GTF
                    (default: stdout)
  -m <min_len>     minimum input transcript length to include in the merge
                    (default: 50)
  -c <min_cov>     minimum input transcript coverage to include in the merge
                    (default: 0)
  -F <min_fpkm>    minimum input transcript FPKM to include in the merge
                    (default: 1.0)
  -T <min_tpm>     minimum input transcript TPM to include in the merge
                    (default: 1.0)
  -f <min_iso>     minimum isoform fraction (default: 0.01)
  -g <gap_len>     gap between transcripts to merge together (default: 250)
  -i               keep merged transcripts with retained introns; by default
                   these are not kept unless there is strong evidence for them
  -l <label>       name prefix for output transcripts (default: MSTRG)

2018-02-13

[][][] SQLからmatplotlibでグラフを作ってFlaskで表示

@app.route('/')
def index():
    import matplotlib.pyplot
    from matplotlib.backends.backend_agg import FigureCanvasAgg
    import io
    import random

    session = Session()
    q = session.query(Diff).filter(Diff.id <6).filter(Diff.id > 1).all()
    fig, ax = matplotlib.pyplot.subplots()
    ax.set_title(u'test')
    for y in q:
        x_ax = ["cont4", "cont5", "cont6", "sample7", "sample8", "sample9"]
        y_ax = [y.Nb_cont04, y.Nb_cont05, y.Nb_cont06, y.Nb_samp07, y.Nb_samp08, y.Nb_samp09]
        ax.plot(x_ax, y_ax, '--o', label = y.id)
        ax.legend(loc='upper right')

    canvas = FigureCanvasAgg(fig)
    buf = io.BytesIO()
    canvas.print_png(buf)
    data = buf.getvalue()

    response = make_response(data)
    response.headers['Content-Type'] = 'image/png'
    response.headers['Content-Length'] = len(data)
    return response

こんな感じで複数の検索項目をグラフ化。

ax.plot(x_ax, y_ax, '--o', label = y.id)

ax.legend(loc='upper right')

この部分で凡例の配置や、項目名、グラフの線やポイントの形などを色々制御できる。

f:id:k-kuro:20180213225346p:image


なお、y軸の数値がintの場合は上のコードでよかったのだが、floatになるとちょっとおかしなことになった。

ValueError: could not convert string to float: "cont4"

というエラーが出てしまう。どうも文字列をFloat型に書き換えようとしてしまってエラーになっている模様。これの解決策として、x軸を一旦なしにして描画(ただしその場合項目が0,1,2...となる)し、あとで文字列へ置き換えてやるコードを追加することになる。

        y_ax = [y.Nb_cont04, y.Nb_cont05, y.Nb_cont06, y.Nb_samp07, y.Nb_samp08, y.Nb_samp09]
        ax.plot(y_ax, '--o', label = y.id)
        matplotlib.pyplot.xticks([0,1,2,3,4,5], ["cont4", "cont5", "cont6", "sample7", "sample8", "sample9"]

こんな感じ。

現状ではグラフを新しいウインドウを開いて表示する形になり、表と見比べるのがやりにくい。表示を同じウィンドウの中に埋め込みたい。

img_data = urllib.parse.quote(data)

return render_template('/rna/exp/index_nb_diff.html', form=form, contents=q, counts=counts, img_data=img_data)

というような感じでurllibでbase64エンコードしてrender_templateでhtmlに引き渡し、templateの方で

<img src="data:image/png:base64,{{ img_data }}"/>

とイメージタグで受け取ると、その場所に埋め込んでもらえるようだ。

とりあえず遺伝子名やGOでSQLに検索をかけてその遺伝子発現パターンを視覚化して評価するアプリが完成した。

次はやはり逆に発現パターンからGOを絞るプログラムが必要だな。ここからはやはりRのおでましかな。

2018-02-12

[][]marplotlibでグラフをWeb上に表示

さて、Pythonでグラフをかけるようになったので、今度はそれをFlaskでweb上に表示させてみたい。

まえに参照したページはpython2系のコードだったようで、そのままではPython3ではエラーが出る。

なので、エラーを一つずつ潰していったところ

@app.route('/')
def index():
    import matplotlib.pyplot
    from matplotlib.backends.backend_agg import FigureCanvasAgg
    import io
    import random

    fig, ax = matplotlib.pyplot.subplots()
    ax.set_title(u'Random')
    x_ax = range(1, 256)
    y_ax = [random.randint(1, 1024) for x in x_ax]
    ax.plot(x_ax, y_ax)

    canvas = FigureCanvasAgg(fig)
    buf = io.BytesIO()
    canvas.print_png(buf)
    data = buf.getvalue()

    response = make_response(data)
    response.headers['Content-Type'] = 'image/png'
    response.headers['Content-Length'] = len(data)
    return response

あとfrom Flask import make_responseもプログラムの先頭あたりで加えておく。

Python2版からの変更点は

cStringIOをimportしていたところをioをimportに変え、

buf = cStringIO.StringIO()

だったところを

buf = io.BytesIO()

に変えた。

これにより、アプリケーションルートにアクセスすれば

f:id:k-kuro:20180212224558p:image

こんな感じに表示されることを確認。

さて、次はSQLのデータをグラフに表示することを目指す。

@app.route('/')
def index():
    import matplotlib.pyplot
    from matplotlib.backends.backend_agg import FigureCanvasAgg
    import io
    import random

    session = Session()
    y = session.query(Diff).filter(Diff.id == 10).one()
    fig, ax = matplotlib.pyplot.subplots()
    ax.set_title(u'test')
    # x_ax = range(1, 256)
    x_ax = ["cont4", "cont5", "cont6", "sample7", "sample8", "sample9"]
    # y_ax = [random.randint(1, 1024) for x in x_ax]
    y_ax = [y.Nb_cont04, y.Nb_cont05, y.Nb_cont06, y.Nb_samp07, y.Nb_samp08, y.Nb_samp09]
    ax.plot(x_ax, y_ax)

    canvas = FigureCanvasAgg(fig)
    buf = io.BytesIO()
    canvas.print_png(buf)
    data = buf.getvalue()

    response = make_response(data)
    response.headers['Content-Type'] = 'image/png'
    response.headers['Content-Length'] = len(data)
    return response

こんな感じでテーブルdiffからクエリを読み込んでグラフ化できた。

f:id:k-kuro:20180212230931p:image

ここまでくれば後はWTFormから検索条件を読んで、それをグラフ化できるようにすれば、なんとなくデータ解析アプリになるな。

2018-02-10

[][][] StringIO

matplotlibでグラフを描画するプログラムの中でcStringIOというモジュールが使われているのだけれど、Python3のモジュールには該当するものが見つからない。

で、調べてみるとPython3では標準モジュールのioにStringIOが含まれているので

from io import StringIO

とするといいらしい。

https://kusanohitoshi.blogspot.jp/2017/01/python3cstringiostringio.html

[][]matplotlib

とりあえずグラフの描写を試してみる。

>>> import matplotlib.pyplot as plt
>>> import random
>>> fig, ax = plt.subplots()
>>> ax.set_title(u'Random')
Text(0.5,1,'Random')
>>> x_ax = range(1, 256)
>>> y_ax = [random.randint(512,1024) for x in x_ax]
>>> ax.plot(x_ax, y_ax)
[<matplotlib.lines.Line2D object at 0x109151fd0>]
>>> plt.show()

python3のコンソールでこのように打つと

f:id:k-kuro:20180211004716p:image

こんな感じにグラフができることを確認した。コレをWeb上に表示させるのが次の課題。

また、複数のグラフを重ねる(複数項目を1グラフに入れる)なら

>>> import matplotlib.pyplot as plt
>>> import random
>>> fig, ax = plt.subplots()
>>> ax.set_title(u'Random')
Text(0.5,1,'Random')
>>> x_ax = range(1, 256)
>>> y_ax = [random.randint(512,1024) for x in x_ax]
>>> ax.plot(x_ax, y_ax)
[<matplotlib.lines.Line2D object at 0x108b890f0>]
>>> x_ax = range(1, 256)
>>> y_ax = [random.randint(128,256) for x in x_ax]
>>> ax.plot(x_ax, y_ax)
[<matplotlib.lines.Line2D object at 0x108bdf898>]
>>> plt.show()

このように続けてプロットするだけでいいらしい。

f:id:k-kuro:20180211010147p:image

2018-02-09

[][][]sqlのDBからグラフを作って表示

集計してテーブルを作って、個々のデータから生データであるBAMを参照することはできるようになったが、テーブルを眺めていてもそれ以上進まないので、テーブルからグラフを描かせていきたい。

https://qiita.com/5t111111/items/3d9efdbcc630daf0e48f

この辺参考になるかな

[]cufflinksのFPKMを遺伝子ごとに集計して、サンプル間で比較する。

常法的にはcuffdiffでFPKMを比較し、cummeRbandでグラフ化するのだけど、サンプル数が多いと処理もけっこう大変になってくる。

とりあえずFPKMを遺伝子ごとに集計して並べて見れるようにしてみることにする。

cufflinksで出力される値に幾つか項目を加えてテーブルを作っている。

class Rnaseq(Base):
    """
    RNA seq database
    """
    __tablename__ = 'rnaseq'
    tracking_id = Column(String(64))
    class_code = Column(String(64))
    nearest_ref_id = Column(String(64))
    gene_id = Column(String(64))
    gene_short_name = Column(String(64))
    tss_id = Column(String(64))
    locus = Column(String(64))
    length = Column(Integer)
    coverage = Column(Integer)
    FPKM = Column(Float)
    FPKM_conf_lo = Column(Float)
    FPKM_conf_hi = Column(Float)
    FPKM_status = Column(String(64))
    id = Column(Integer, primary_key=True)
    sample_id = Column(Integer, ForeignKey('samples.id'))
    xtr_id = Column(Integer, ForeignKey('rnaseq_xtr.id'))
    gene_id2 = Column(Integer, ForeignKey('Genes.id'))

gene_id2でjoinするGenesテーブルの方は

class Genes(Base):
    __tablename__ = 'genes'
    id = Column(Integer, primary_key=True)
    a_id = Column(String(64))
    b_id = Column(String(64))
    gene_id = Column(String(64))
    rnaseq = relationship('Rnaseq', backref='genes', lazy='dynamic')

そこにさらに比較用テーブルを作成して

class Diff(Base):
    __tablename__ = 'diff'
    id = Column(Integer, primary_key=True)
    gene_id = Column(String(64))
    samp1 =  Column(Float)
    samp2 =  Column(Float)
    samp3 =  Column(Float)
    ......
    ......
    samp1l =  Column(String(64))
    samp2l =  Column(String(64))    
    samp3l =  Column(String(64))
    ......
    ......

こんな感じに用意してgene_idにgenesテーブルのb_id、samp1〜に各サンプルのFPKM、samp1l〜にlocusを割り当ててコピーする。

コピーの仕方は

q = select(["genes.id, genes.b_id, rnaseq.FPKM, rnaseq.locus"],
                and_("genes.b_id = rnaseq.gene_id", "rnaseq.sample_id = 1"),
                from_obj=["genes, rnaseq"]) 
ins = insert(Diff).from_select(
                (Diff.id,
                Diff.gene_id,
                Diff.samp1, 
                Diff.samp1l),
                q)
conn.execute(ins)
q = session.query(Genes.id, Genes.b_id).filter(not_(Genes.b_id.in_(session.query(Rnaseq.gene_id).filter(Rnaseq.sample_id = 1)))
ins = insert(Diff).from_select(
                (Diff.id,
                Diff.gene_id),
                q)
conn.execute(ins)

まずsample1をコピー。つぎにこれにsample2を加えてdiffテーブルのコピーであるdiff_tempにinsert

つぎにsample3を加えてdiffにinsertしなおし・・・と順々に行ったり来たりを繰り返す。ただし、このままでは各sampleに含まれているgeneだけに絞られていってしまうので、sampleに含まれていなかったぶんはその都度前のテーブルからそのままinsertで補充していく必要がある。

2018-02-08

[][]Nbのリファレンス

$ 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になってしまっている。なんでかなー。

f:id:k-kuro:20180208093923p:image

アプリは突貫で一応できており、検索、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の4つが乗っているが、変換したgtfではexonしか残っていないな。これがマッピングでアノテーションが全然反映されていない原因かな。

f:id:k-kuro:20180208100552p:image

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の純正ツールなんだから間違いないだろう。

さて、状況はすべて把握できた。いざマッピング。

f:id:k-kuro:20180208130018p:image

バッチリじゃなかろうか。とりあえず出来上がった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ファイルを作成し、リファレンスに指定してやること。

2018-02-07

[][] 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カラム。コレはいかに?

[][]リファレンスとアノテーションファイルの中身を理解する。

とりあえずちゃんとできているArabidopsisで使っているファイルを見てみると

  • マッピングに使うリファレンス:Arabidopsis_thaliana.TAIR10.dna.toplevel.fa

(IGVのreferenceにもそのまま使う)

$ head Arabidopsis_thaliana.TAIR10.dna.toplevel.fa
>1 dna:chromosome chromosome:TAIR10:1:1:30427671:1 REF
CCCTAAACCCTAAACCCTAAACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAAT
CTTTAAATCCTACATCCATGAATCCCTAAATACCTAATTCCCTAAACCCGAAACCGGTTT
CTCTGGTTGAAAATCATTGTGTATATAATGATAATTTTATCGTTTTTATGTAATTGCTTA
TTGTTGTGTGTAGATTTTTTAAAAATATCATTTGAGGTCAATACAAATCCTATTTCTTGT
GGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAA
GCTTTGCTACGATCTACATTTGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTT
ATCTCAAGAATCTTATTAATTGTTTGGACTGTTTATGTTTGGACATTTATTGTCATTCTT
ACTCCTTTGTGGAAATGTTTGTTCTATCAATTTATCTTTTGTGGGAAAATTATTTAGTTG
TAGGGATGAAGTCTTTCTTCGTTGTTGTTACGCTTGTCATCTCATCTCTCAATGATATGG

つまりクロモソームごとのゲノム配列そのもの。

  • cufflinks, cuffdiffで使うアノテーションファイル:Arabidopsis_thaliana.TAIR10.38.gtf
$ head Arabidopsis_thaliana.TAIR10.38.gtf
#!genome-build TAIR10
#!genome-version TAIR10
#!genome-date 2010-09
#!genome-build-accession GCA_000001735.1
#!genebuild-last-updated 2010-09
1	araport11	gene	3631	5899	.	+	.	gene_id "AT1G01010"; gene_name "NAC001"; gene_source "araport11"; gene_biotype "protein_coding";
1	araport11	transcript	3631	5899	.	+	.	gene_id "AT1G01010"; transcript_id "AT1G01010.1"; gene_name "NAC001"; gene_source "araport11"; gene_biotype "protein_coding"; transcript_source "araport11"; transcript_biotype "protein_coding";
1	araport11	exon	3631	3913	.	+	.	gene_id "AT1G01010"; transcript_id "AT1G01010.1"; exon_number "1"; gene_name "NAC001"; gene_source "araport11"; gene_biotype "protein_coding"; transcript_source "araport11"; transcript_biotype "protein_coding"; exon_id "AT1G01010.1.exon1";
1	araport11	CDS	3760	3913	.	+	0	gene_id "AT1G01010"; transcript_id "AT1G01010.1"; exon_number "1"; gene_name "NAC001"; gene_source "araport11"; gene_biotype "protein_coding"; transcript_source "araport11"; transcript_biotype "protein_coding"; protein_id "AT1G01010.1"; protein_version "1";
1	araport11	start_codon	3760	3762	.	+	0	gene_id "AT1G01010"; transcript_id "AT1G01010.1"; exon_number "1"; gene_name "NAC001"; gene_source "araport11"; gene_biotype "protein_coding"; transcript_source "araport11"; transcript_biotype "protein_coding";

左から順に

1 ->chromosome

araport11 ->データソース

gene -> feature

3631 ->start pos

5899 ->end pos

. ->score

  1. ->direction

. ->frame (0,1,2 or .)

gene_id "AT1G01010"; gene_name.... -> attribute

gffファイルだとattributeの書き方が違う

  • igv.jsなどでreference配列の表示に使う:Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai
$ head Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.fai
1	30427671	55	60	61
2	19698289	30934909	60	61
3	23459830	50961558	60	61
4	18585056	74812441	60	61
5	26975502	93707303	60	61
Mt	366924	121132452	60	61
Pt	154478	121505547	60	61

クロモソームごとの indexが作られている。カラムは左から

NAME Name of this reference sequence

LENGTH Total length of this reference sequence, in bases

OFFSET Offset within the FASTA file of this sequence's first base

LINEBASES The number of bases on each line <- 行の折り返し

LINEWIDTH The number of bytes in each line, including the newline

OFFSETとは実際の配列の長さに加え、ヘッダにある文字や、改行をカウントしていってファイルの先頭から何バイトがそのseqの先頭であるかを示すもの。

  • igv.jsでアノテーションのトラックを表示するファイル: TAIR10.bed
$ head TAIR10.bed
chr1	3630	5899	AT1G01010.1	0	+	3759	5630	0	6	283,281,120,390,153,461,	0,365,855,1075,1543,1808,	NAC001|ANAC001	NAC domain containing protein 1
chr1	5927	8737	AT1G01020.1	0	-	6914	8666	0	10	336,633,76,67,86,74,46,90,48,167,	0,509,1229,1456,1636,1834,2014,2308,2489,2643,	ARV1	Arv1-like protein
chr1	6789	8737	AT1G01020.2	0	-	7314	8666	0	8	280,294,86,74,46,90,48,167,	0,367,774,972,1152,1446,1627,1781,	ARV1	Arv1-like protein
chr1	11648	13714	AT1G01030.1	0	-	11863	12940	0	2	1525,380,	0,1686,	NGA3	AP2/B3-like transcriptional factor family protein
chr1	23145	31227	AT1G01040.1	0	+	23518	31079	0	20	1306,114,211,395,220,173,123,161,234,151,183,162,96,629,98,191,906,165,407,326,	0,1396,1606,1895,2378,2679,2935,3146,3397,3716,3953,4226,4472,4657,5562,5744,6014,7001,7264,7756,	EMB76|SIN1|SUS1|ATDCL1|DCL1|ASU1|EMB60|CAF	dicer-like 1
chr1	23415	31120	AT1G01040.2	0	+	23518	31079	0	20	1036,114,211,395,220,173,123,161,234,151,183,165,96,629,98,191,906,165,407,219,	0,1126,1336,1625,2108,2409,2665,2876,3127,3446,3683,3956,4202,4387,5292,5474,5744,6731,6994,7486,	EMB76|SIN1|SUS1|ATDCL1|DCL1|ASU1|EMB60|CAF	dicer-like 1
chr1	28499	28706	AT1G01046.1	0	+	28499	28499	0	1	207,	0,	MIR838A	MIR838a; miRNA
chr1	31169	33153	AT1G01050.1	0	-	31381	32670	0	9	255,82,121,66,108,66,29,124,125,	0,351,523,763,918,1112,1261,1377,1859,	AtPPa1|PPa1	pyrophosphorylase 1
chr1	33378	37757	AT1G01060.3	0	-	33991	37061	0	10	211,347,1074,81,234,62,112,181,26,189,	0,602,1022,2188,2351,3245,3431,3644,3994,4190,	LHY1|LHY	Homeodomain-like superfamily protein
chr1	33665	37780	AT1G01060.2	0	-	33991	37061	0	8	662,1074,81,234,62,112,181,408,	0,735,1901,2064,2958,3144,3357,3707,	LHY1|LHY	Homeodomain-like superfamily protein

左から

chr1 ->クロモソーム番号

3630 ->start pos

5899 ->end pos

ここまで必須

AT1G01010.1 ->名前

0 ->score

'+ ->方向

3759 ->CDSのstart pos

5630 ->CDSのend pos

0 ->なにもなし

6 ->exonの数

283,281,120,390,153,461, ->各exonのサイズ

0,365,855,1075,1543,1808, ->各exonの先頭pos

ここまでオプション

LHY1|LHY

Homeodomain-like superfamily protein

この2つは独自拡張しているのかもしれない。

方向までの6カラムで一応用をなすのでそれ以外はigvにはあってもなくてもいいというわけか。なお、カラムの並びが違うだけでgtf/gffも同じ内容を持っているので、igvでも指定さえすればどれでも使えるということだな。

2018-02-02

[][] arabidopsisでhisat2

さて、これまでヒトゲノム研究でRNAseqパイプラインをせっせと作っていたのだが、縁あって植物科学の世界に舞い戻ってきた。

そこで植物のRNAseqを解析することになったのだが、やはりと言うか、色々と環境が整備されていなくて、いちいち対応していく必要がある。

前任者、というかコラボでやってもらっていたデータをはいっと渡されはしたものの、どういうストラテジーで解析されたファイルなのかわからないまま眺めてみても、やはり今ひとつ状況が把握できない。

ということで、まずはシングルエンドのシーケンスリードをtrimmomaticで処理して、Hisat2にかけてマッピングまでやってみることにする。

trimmomaticのシングルエンドリード処理はやったことがなかったのでパイプラインを見直して、まずは対応させるところから。

基本的にはSEというオプションを付けて、入力するファイル、出力するファイル名を1つづつにしてやればそれで良いようだ。

次にHisat2でのマッピング。

そもそもリファレンスファイルから準備しなくてはならない。

http://plants.ensembl.org/info/website/ftp/index.html

こちらのサイトに植物のゲノム情報が集約されているようなので、そこから該当するファイルを貰ってくることにした。

とりあえずはarabiのDNAとアノテーションgtfファイルをダウンロードする。

更にHisatでマッピングするにはインデックスファイルが必要となるのだが、Hisat2のページで配布されているのは動物関係しかないようなので、これも自前で作成してやらねばならない。

ダウンロードしたfaファイルの階層にて

$ hisat2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Arabidopsis_thaliana.TAIR10.dna.toplevel

こういうふうに打つと、インデックスファイルが8分割して生成した。

Headers:
    len: 119481543
    gbwtLen: 119481544
    nodes: 119481544
    sz: 29870386
    gbwtSz: 29870387
    lineRate: 6
    offRate: 4
    offMask: 0xfffffff0
    ftabChars: 10
    eftabLen: 0
    eftabSz: 0
    ftabLen: 1048577
    ftabSz: 4194308
    offsLen: 7467597
    offsSz: 29870388
    lineSz: 64
    sideSz: 64
    sideGbwtSz: 48
    sideGbwtLen: 192
    numSides: 622300
    numLines: 622300
    gbwtTotLen: 39827200
    gbwtTotSz: 39827200
    reverse: 0
    linearFM: Yes
Total time for call to driver() for forward index: 00:02:07

まあこんなもんでしょう。

さてさて、次にHisat2にかけてみよう。

そもそもfastqファイルが4分割されているようなんで、これをどの段階でマージすればいいのかがわからず。

そのまま以前作ったパイプラインに流してみると、どうもBAMファイルに変換して、ソートをかけようとするとエラーが出て先に進まないようなので、ソートする前、もしくはBAMに変換する前にどうにかしないといけないのかもしれない。


やってみたところ、BAMに変換後、samtools mergeで結合したBAMファイルに対してsort-indexとかければちゃんと動作し、出来上がったBAMファイルをIGVブラウザで閲覧することができた。

f:id:k-kuro:20180202162559p:image

つぎに、これをcufflinksにかけて発現量を計算してみる。

$ cufflinks -p 4 -q -g "$gtf" "$output_dir"/"$exp".bam -o "$output_dir2"

gtfファイルやインプット・アウトプットは前もって指定してあるのでこんな感じ。4コアなiMacでうりゃっとかけると

f:id:k-kuro:20180202164242p:image

見事に4コア使い切っていて清々しい。1サンプル40分ってところだね。(2.7 GHz Intel Core i5 4core)