不求东西

“不求南北,各奔东西”

基于SQLAlchemy构建精简的API

SQLAlchemy是一个强大的数据库操作工具,也是目前Python的ORM实际上的标准。它的便利性主要体现在以下三点:

  1. 它支持大多数的主流数据库,如SQLite,Oracle,MySQL,Microsoft SQL server等,为不同的数据库提供了统一的接口。
  2. 它提供了一种面向对象的基于SQL语言的抽象(也就是ORM),这使得我们操作数据库时候,可以写出更干净,也更Pythonic的代码。免除了SQL语言的学习成本,也避免了将SQL代码和Python代码混合起来写所带来的复杂度。
  3. 对于底层性能的优化(如数据库连接池的创建),SQLAlchemy都会自动完成。

问题的由来

虽然SQLAlchemy很好很强大,可是有时候我也会感到它的繁琐之处。比如要完成一项查询,查找Dataset对象中的两列:id和name,过滤条件是id列的值小于100:

results=DBSession.query(Datasets).filter(Datasets.id<100).all()
map(lambda x:[x.id,x.name],results)

虽然只有两行,可是冗余的部分还是有的,我尝试将不变的部分剥离出来,构建了一个精简化的函数:

def meta_API(ORM_entity,selector,filterer):
    if isinstance (ORM_entity,list):
        return map(selector,
                   filterer(DBSession.query(*ORM_entity)))
    else:
        return map(selector,
                   filterer(DBSession.query(ORM_entity)))

这样就把selector(选择器)和filterer(过滤器)这两个过程抽象了出来,于是可以将前面对Datasets对象的查询过程改成:

meta_API(ORM_entity=Datasets,
         selector=lambda x:[x.id,x.name],
         filterer=lambda x:x.filter(Datasets.id<100).all())

到底是看上去更清晰了,还是更古怪了,这倒是一个见仁见智的问题。

构建更高层的API

以meta_API为接口,可以写出更多的语法糖,比如dump(输出全部),head(输出首行),search(按表达式搜索)以及joint(连接表)等,如下:

def meta_dump_API(ORM_entity,selector,filterer=lambda x:x):
    """
    Return all the result after filtering and selecting.

    The function is a syntactic sugar derived from meta_API. You can skip `.all()` in the filterer.
    """
    return meta_API(ORM_entity=ORM_entity,
                    filterer=lambda x:filterer(x).all(),
                    selector=selector)

def meta_search_API(ORM_entity,selector,predictives):
    """
    Return all the result matches the predictives after selecting.

    The function is a syntactic sugar derived from meta_API. You can use predictives instead of filterer to skip `.filter().all()`.
    """
    if isinstance (ORM_entity,list):
        return meta_API(ORM_entity=ORM_entity,
                        filterer=lambda x:x.filter(*predictives).all(),
                        selector=selector)
    else:
        return meta_API(ORM_entity=ORM_entity,
                        filterer=lambda x:x.filter(predictives).all(),
                        selector=selector)

def meta_head_API(ORM_entity,selector,filterer=lambda x:x):
    """
    Return the first result after filtering and selecting.

    The function is a syntactic sugar derived from meta_API. You can skip `.one()` in the filterer.
    """
    return meta_API(ORM_entity=ORM_entity,
                    selector=selector,
                    filterer=lambda x:filterer(x).one())

def meta_joint_API(ORM_entities, filterer, selector,
                   inner=Dataset,outter=None):
    """
    Return the result after joining,filtering and selecting.

    The function is a syntactic sugar derived from meta_API. You can skip `.join()` in the filterer.
    """
    return meta_API(ORM_entity=[inner]+ORM_entities,
                    selector=lambda x:selector(x[1:]),
                    filterer=lambda x:filterer(x.join(outter)))

可以看出,meta_API这一抽象还是有其实用性的。

定制API

最后就可以根据特定数据库的数据结构自己定制API了,以我们实验室的内部数据库为例:

from sqlalchemy import engine_from_config
from settings import DATABASE_CONFIG
engine=""
def initDB():
    global engine;
    engine=engine_from_config(DATABASE_CONFIG,"")
    DBSession.configure(bind=engine)
def dump_CellLine(filterer=lambda x:x):
    return meta_dump_API(ORM_entity=CellLine,
                         selector=lambda x:(x.id,x.name),
                         filterer=filterer, )

def head_CellLine():
    return meta_head_API(CellLine,
                         lambda x:(x.id,x.name))
def joint_CellLine(selector=lambda):
    return meta_joint_API(ORM_entities=[CellLine,Dataset],
                          selector=lambda x:x,
                          filterer=lambda x:x.all(),
                          inner=Dataset,
                          outter=CellLine)

注:Dataset中需要加ForeignKey参数

class Dataset(Base):
    __tablename__='mytable_datasets'
    id = Column(Integer, primary_key=True)
    cell_line_id = Column(Integer,ForeignKey("mytable_celllines.id"))

感觉joint的那个API设计得还是很繁琐,不过我目前只能想到这么多。

中文化的Galaxy

国内应该有相当一部分人用过Galaxy吧,我不久前给Galaxy项目提交了一份中文化翻译的模版,可将官方Galaxy的网页界面部分汉化,显示效果如图:


截图中的字体不一致的现象是我的浏览器不给力造成的。。

该汉化版本已被接收,想使用它的话,进入官方Galaxy页面,然后将你所用浏览器的语言偏好设为中文,刷新页面即可看到效果。
目前中文化的完成度不高,希望在一年后能看到更好的中文Galaxy界面。

注:Galaxy是专为生物学家设计的在线数据分析平台,开发人员来自于psu和Emory University,它的对手是Broad Institude开发的GenePattern

使用Aspera从EBI或NCBI下载基因组数据

做基因组数据分析,可能经常从NCBI的GEO/SRA或者EBI的ENA数据库下载高通量的数据,动辄几十G的数据用wget下载实在太纠结,这时就要用到神器-Aspera了。

使用Aspera,最简单的方法当然就是使用浏览器插件Aspera Connect了,跟迅雷、Flashget的用法差不多,直接单击Aspera支持的下载地址,就自动切换到Aspera的窗口开始下载了。

当我们登录到自己的服务器终端里面的时候,可能更希望在终端里直接下载数据,而不是先把数据下载到自己的硬盘里,再上传到服务器,这种情况下带有窗口界面的Aspera Connect就无法使用了吗?

当然可以,Aspera Connect安装包里内置了Aspera的命令行工具,这里对其安装和使用方法简要介绍一下:

安装

首先,到aspera网站下载你的操作系统对应的aspera connect。(如果选Linux,下载以后会是一个几M大,内嵌二进制代码的shell脚本。。) 。不需要root或者sudo权限,直接安装之:
$ sh aspera-connect-2.4.7.37118-linux-64.sh
安装好以后,会在HOME目录下新建一个叫.aspera的目录,有两个文件比较重要:
一个是ascp的可执行文件:
~/.aspera/connect/bin/ascp
另一个ascp的密钥文件:
~/.aspera/connect/etc/asperaweb_id_dsa.putty

建议将密钥备份到HOME目录下方便使用:
$ cp ~/.aspera/connect/etc/asperaweb_id_dsa.putty ~/
再把aspera-license复制到系统目录
~/.aspera/connect/etc$ sudo cp aspera-license /usr/local/bin/
再把ascp可执行文件的路径加入PATH变量中,或者将其拷贝到当前目录。

使用

执行以下两条命令(注意最后要加点号“.”,表示当前目录)
从EBI下载:
$ ascp -i ~/asperaweb_id_dsa.putty era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA012/ERA012008/sff/library08_GJ6U61T06.sff .
从NCBI下载:
$ ascp -i ~/asperaweb_id_dsa.putty anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/litesra/SRR/SRR096/SRR096072/SRR096072.lite.sra .

这个时候的速度相比于wget,应该已经很快了,大约能达到9Mb/s以上,如果还嫌慢,可以在-i 参数的前面添加几项设置,像这样:
ascp -QT -l 100M -i ~/asperaweb_id_dsa.putty era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA012/ERA012008/sff/library08_GJ6U61T06.sff .这样可以将速度提高到20Mb/s左右,偶尔能达到100Mb/s。

ascp下载地址的获取

以EBI上的SRR346368这套数据为例。首先到EBI页面里,找到你想要下载的文件,将指针移到这个文件的”ftp”这一列,即可看到其ftp地址,例如: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR346/SRR346368/SRR346368.fastq.gz,
然后呢:将 ftp://ftp.sra.ebi.ac.uk 换成 era-fasp@fasp.sra.ebi.ac.uk即可:
$ ascp -i ~/asperaweb_id_dsa.putty era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR346/SRR346368/SRR346368.fastq.gz .
NCBI的SRA数据库也是同样的方法,即可获取其ascp下载地址。

小技巧

如果嫌每次都输入密码太麻烦,可以在命令行或.profile中设置ASPERA_SCP_PASS这个环境变量:
export ASPERA_SCP_PASS=你的aspera密码
即可。

更多的说明

请参见官方的SRA下载手册:
NCBI: http://www.ncbi.nlm.nih.gov/books/NBK47540/
EBI: http://www.ebi.ac.uk/ena/about/sra_data_download

P.S. 好像我身边的人喜欢EBI的偏多,因为对于同一套公开数据,EBI可以下载到原始测序文件(Fastq格式),而NCBI中只能下载压缩转换格式后的文件(SRA格式)。不知道大家喜欢用哪一种?

一一年读书小记(技术篇)

Python学习手册 第四版 -Mark Lutz 著

  • 直白,不绕弯子,同时又耐心地突出重点,是我喜欢这本书的原因。
  • 这本书的前20章我认为是目前Python入门书里面最好的,没有之一。
  • 而本书的后半部分(类和OOP,高级话题),不知道是因为作者难于把握,还是因为译者水平太洼,和前面的内容相比,总觉得欠点火候。

Python Cookbook 第二版 -Alex Marteli 编

  • 正如其名,这本书像菜谱一样,可以从任意一页开始看,很多小节的标题看上去很有趣,而且讲得也够深入。
  • 当初希望把这本书当成Python文集来看,却发现不如想象中那么好:不少内容对我来说都是屠龙之技,日常生活基本用不到,比如说Tkinter界面编程、时间计算、分布式编程等等。
  • 不过当遇到相关问题时,翻一翻这本书还是有帮助的。

黑客与画家 -Paul Graham 著

  • 话说这本书2004年就被O’Reilly出版了,可是中文版到了2011年才出,可见国内出版业的滞后性。。
  • 看看作者的学历:本科在康奈尔修哲学,博士期间在哈佛搞人工智能,读完了博士又去欧洲学绘画,后来再后来呢,就发达了(不知道的Google下)。总之确实让人感慨:要想成点事,这人还真得能折腾。
  • 这书写的确实挺犀利的,至少给hacker们,或者给那些向往hacker的程序员们,提供了一个很好的心理优越感的来源。
  • 推荐一看,再不济,也可以给那些没接触过Lisp的人扫扫盲,何况这本书从头到尾都挺耐读的。

人月神话 -Frederick P. Brooks 著

  • 这本书的观点在当时也算是一针见血吧,让我印象深刻的观点有:越帮越忙的“焦油坑”,每个人或多或少幻想过的“银弹技术”,以及编程过程中任何“高级工具”和“编程方法”无法消除的“本质的复杂度”等。
  • 而对于那些曾经在项目中有挫折感的程序员来说,没有比这本书更好的励志读物了。

人件 – Tom DeMarco/Timothy Lister 著

  • 看完这本书以后,我立刻换了座位,从整天充满闲聊鬼扯而且拥挤不堪的过道上,搬到了一个相对宽敞,相对安静的角落里。
  • 读后感:对于个人来说,环境舒适度(噪音、面积)很重要。对于团队来说,内部精英意识的形成很重要

Shell脚本学习指南 第二版 – Arnold Robbins/Nelson H.F.Beebe 著

  • 写得很细致,介于使用手册与教程之间,把bash中所有命令和参数都写了一遍,但说老实话,有点枯燥。
  • 如果不是系统管理员,能用到其中讲的10%就非常不错了。

实用Common Lisp编程 -Peter Seibel著

  • 真心感觉这本书没有介绍的那么好,难道是翻译的问题?
  • 翻遍整本书,也没找到这个很基本的问题的答案:怎样能不在REPL环境中,而是在命令行下直接运行lisp脚本。对于lisp解释器在命令行下的用法方面的介绍,作者好像惜字如金。对于一本以“实用”作为噱头的书来讲,实在是令人失望。
  • 本来想用Lisp代替Python和awk的,这本书看了六章以后还是放弃

Refseq与Gene Symbol之间的区别

Refseq

  • 目的:为中心法则中自然存在的分子,提供唯一的命名。
  • 格式:目前做基因组学最最常见的是NM_XXX(mRNA),NP_XXX(蛋白质),NC_XXX(染色体/基因组/序列)

Gene Symbol

  • 目的:为每一个基因,提供唯一的命名(一般是缩写形式)。
  • 格式:大写字母,后面可加数字

 

区别很明显了,Refseq为分子命名,Gene Symbol为基因命名。

不过有时候这种界限也不是很明显,有时候也可以用mRNA的Refseq名称来对应其基因。

所以有些基因注释文件用Gene Symbol来表示基因,而用Refseq来表示该基因的转录本中不同的mRNA分子(isoform)。

在这种情况下,一个Gene Symbol会对应多个Refseq分子。

ISDSB小记—第三天

今天会议的主题是计算生物学,小记下今天的2个演讲人吧。

Yijun Ruan:ChIA-PET技术大牛

  • ChIA-PET,说成是ChIP-Seq技术的升级版大概也不为过吧
  • ChIA-PET不仅可以观察那些区域是转录因子的结合位点,还可以推测出到这些区域附近的三维结构,后一点是ChIP-Seq目前做不到的
  • Yijun Ruan他们做了ERα,CTCF和PolII三个蛋白质的ChIA-PET

Yingrui Li (李英睿):BGI新星

  • 穿着:短裤+拖鞋
  • 谈到了个人测序,不光要测SNP
  • 北大04级生物系,大二入华大,本科毕业前以第一作者发Nature文章,现为BGI某团队负责人

看过这哥们的简历,我是觉得,现在每天吭哧吭哧地干活,幻想着七八年后能把文章发到Nature、Science上,赶上大二时候的他,这样的想法确实挺二的。

ISDSB小记—第一天

今天大概有二十个报告人,个人水平有限,小记几个熟悉的吧。

Michael Levine:果蝇胚胎细胞中“影增强子”对转录机制的调节。

  • “影增强子” (shadow enhancer)是相对“主增强子” (primary enhancer)而言的:“主增强子”距离靶基因 (target gene)较近,而“影增强子”则较远
  • “影增强子”可能在发育过程中,保证了基因表达的精确性和可重现性 (reproducibility)
  • “影增强子”可能与Pol II 的停顿 (pause)和释放 (release)有关

Micheal Q. Zhang:用于ChIP-seq数据的贝叶斯模型。

  • 用于从ChIP-seq实验的结果中,获取显著的转录因子 (transcription factor)的结合区域,这一过程也称作peak calling
  • 与熟悉的MACS软件进行了比较

Shirley Liu:染色质动力学中转录因子结合的独特模式

  • 在雄性激素 (Androgen receptor)的刺激下,FoxA1转录因子结合在特定的DNA序列上
  • 在此过程中,结合位点处的转录因子会把附近的核小体挤走

P.S. 大早上就和几个同学去美兰湖会场蹭会议。可下午3点钟不到,认识的同学就走得差不多了。会议终于在晚上7点钟结束了,走出会场却发现,外面电闪雷鸣地下着大雨,无奈,冒雨乘地铁独归。

sed, Awk与Perl单行脚本快速入门

今天偶然看到了这几篇文章,虽然不是很懂,但还是觉得将来会用到,分享一下:

sed的常用单行脚本:

http://sed.sourceforge.net/sed1line_zh-CN.html

Awk的等价实现:

http://linuxtoy.org/archives/sed-awk.html

Perl的等价实现:

http://marlonyao.iteye.com/blog/1013248

这些很简单的功能,用Python完成起来可就稍微复杂了点吧,虽然比可读性的话其实还是Python (废话~)。

不过我倒觉得Perl有些设计确实很巧妙 (比如-pe和-ne这些命令行选项,很适合写单行脚本),有空还真想鼓捣一下Perl。