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で補充していく必要がある。