搜索
查看: 2768|回复: 1

[Other] 生信编程直播第10题的 SQL 解决方案

[复制链接]

2

主题

3

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-2-21 14:28:52 | 显示全部楼层 |阅读模式
原题

生信编程直播第10题:根据指定染色体及坐标得到位置信息
这个题目是来源于生信编程直播学员的真实需求! 当然,我也遇到过这个需求:http://www.bio-info-trainee.com/1001.html
任意给定基因组的 chr:pos, 判断它在哪个基因上面?这个程序难吗?
基因的chr,start,end都是已知的 (这个文件需要下载,可以是UCSC,或者NCBI)
学术一点讲述这个问题:已知CNV数据在染色体上的position如chr1:2075000-2930999,怎样批量获取其对应的Gene Symbol呢(批量)
而且还可以更复杂,我不仅想知道自己指定的坐标是否在某个基因上: 如果在基因上,我还想知道在内含子还是外显子,进一步说是第几个外显子或者内含子! 如果不在基因上,我想知道我指定的坐标距离哪一个基因的TSS(转录起始位点)最近呢?
就是对指定位点或者区间做genomic features的注释!
当然,这个大把的已经造好的轮子可以做到,比如annovar,snpEFF,VEP就是做单个位点注释genomic features的,其它CHIP-seq的软件,都可以做单个区间的批量注释到genomic features!
既然有现成的轮子,我推荐大家去学习一下就好了。
我们的练习题还是简单一点好!
就是指定位置注释是否在基因上面吧!

使用 SQL 语句注释基因

我这里给出的解题方案,是基于 SQL 语句的。

实际上以前我只看 《Perl入门》小骆驼时,也不知道有其他方案,题主用过的这个办法(http://www.bio-info-trainee.com/1001.html)我也曾经用过,原因是这种方法是小骆驼一书贫乏的知识点中(数组列表、哈希表、正则表达式),和这个问题最接近的一种解决方案。
但这个方案缺点同样十分明显,数据量千万不能太大。于是后来我发现了 bedtools,就用起了这个轮子。

后来我又发现,其实如果拿到的不是基因区间,而是地理位置、经纬度的大规模注释,比如手机打车寻找最近出租车的这个问题,基于哈希就没办法算了——总不能把地球所有地方都放进内存里面吧?他们怎么做?总不能也用 bedtools 吧? 李开复 2006 年曾经写过一篇文章《算法的力量》,明确提到过  rtree结构适用于地理位置定位,那么我们能否用 rtree 去做一下试试呢?

有问题问百度。我们如果去百度/Google一下rtree解决实际问题的源码,会发现实际上代码很少。这么重要的问题怎么可能搜不出?原因是他太常见了。和 rtree 同属常见数据结构的 hash 表不同的是,hash 表代码也不好找,但他使用起来太常见了, perl python是直接支持,C++ 有标准库和 Boost 库支持,不需要自己写。rtree虽然 C++ perl python 这些不直接支持,但是实际上他出现在类似这种地方:
SQL 语言的表
R 语言的 dataframe
Python Pandas 包计算包系列

所以,以后拿到类似区间文件去找 overlap,最好的办法,就是把数据放进数据库中(内存指针+存在硬盘里)或者直接通过Python/R 读进内存、以 DataFrame 形式调用。


这里以 SQL 举例,说明如何用 SQL 解决这个问题:

首先

[Python] 纯文本查看 复制代码
import MySQLdb



这里可能会有错误说 import error,没有这个包。如果你安装的是 anaconda 版本的 Python,请在下面多打一行命令,让 notebook 在后台安装:

[Python] 纯文本查看 复制代码
!conda install mysql-python


仍然有问题的话,可能就需要去 yum linux 源库(以CentOS Redhat系统为例)。可以参考我用 linux 裸机部署生物信息分析(RNA/ChIP-Seq)的代码。

成功 import 后,我们连接服务器。通过 Python 的接口连接 UCSC 服务器:


[Python] 纯文本查看 复制代码
db = MySQLdb.connect(host="genome-mysql.cse.ucsc.edu", db='hg38', user="genome")
c = db.cursor()


有人会问为什么连接这里?
首先,UCSC 的 host 和 用户 的相关说明网页在这里: http://genome.ucsc.edu/goldenPath/help/mysql.html
其次,我们在 UCSC 下载基因注释文件,如最新的人类的注释: http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/
会看见有两个 refGene 文件:
refGene.sqlrefGene.txt.gz
我们通常下载的都是第二个,下载完成,解压缩,用 perl 逐行读取。至于第一个,也不知道是干嘛的。
实际上第一个是基因注释文件表的表头。文件主要内容如下:

[SQL] 纯文本查看 复制代码
CREATE TABLE `refGene` (
  `bin` smallint(5) unsigned NOT NULL,
  `name` varchar(255) NOT NULL,
  `chrom` varchar(255) NOT NULL,
  `strand` char(1) NOT NULL,
  `txStart` int(10) unsigned NOT NULL,
  `txEnd` int(10) unsigned NOT NULL,
  `cdsStart` int(10) unsigned NOT NULL,
  `cdsEnd` int(10) unsigned NOT NULL,
  `exonCount` int(10) unsigned NOT NULL,
  `exonStarts` longblob NOT NULL,
  `exonEnds` longblob NOT NULL,
  `score` int(11) DEFAULT NULL,
  `name2` varchar(255) NOT NULL,
  `cdsStartStat` enum('none','unk','incmpl','cmpl') NOT NULL,
  `cdsEndStat` enum('none','unk','incmpl','cmpl') NOT NULL,
  `exonFrames` longblob NOT NULL,
  KEY `chrom` (`chrom`,`bin`),
  KEY `name` (`name`),
  KEY `name2` (`name2`)
) ENGINE=MyISAM DEFAULT CHARSET=latin1;

这里的主要信息告诉我们,这个表的名称是 refGene,然后他有若干列,以及每一列的数据类型是什么。
特别的,name2 这列写的是基因的名称。
OK,那我们写一句 SQL,找位于某一染色体上、 txStart 与 txEnd 之间的所有基因。
[Python] 纯文本查看 复制代码
sql_query = """
SELECT DISTINCT(name2) FROM refGene
     WHERE ((chrom = "%s") AND (txStart >=%d) AND (txEnd < %d)) LIMIT 4
""" 
chrom, beg, end = ["chrX", 44873173, 44874173]
c.execute(sql_query % (chrom, beg, end))
while 1:
    result = c.fetchone()
    if result != None:
        print(result)
    else:
        break


('KDM6A',)

结果被调了出来。即 KDM6A 在 hg38 上位于 chrX:44873173-44874173 这一区间当然,如果这个表过大,需要密集查询,实际上是有本地化的必要的,之前介绍的 pySparkSQL 就是一个威力巨大的工具。但 SparkSQL 引擎装起来还是太麻烦,想简单,也可以试试 sqlite 。






上一篇:rMATS 报错:AttributeError: 'csamtools.AlignedRead' object has no attribut...
下一篇:R语言可视化学习
回复

使用道具 举报

5

主题

37

帖子

485

积分

中级会员

Rank: 3Rank: 3

积分
485
发表于 2017-2-28 15:49:11 | 显示全部楼层
bedtools注释CNV区间基因,挺方便的
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 00:09 , Processed in 0.041606 second(s), 31 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.