搜索
查看: 5772|回复: 4

从Fasta文件取序列的深度解析

[复制链接]

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-2-28 01:57:45 | 显示全部楼层 |阅读模式
本帖最后由 tsznxx 于 2017-2-28 01:57 编辑

本人没有太多的时间,只是偶尔追下大神们的生信课,有时看到些课程,觉得有些可以跟大家分享的就单独拿出来说一下。

详细讲解如何抽取fasta文件里面指定序列名的序列说开去。从Fasta文件里面取序列是新手练script的基本课程,通常的方式用是grep的方式找到对应的序列名称,然后读后面的序列,取出相应的片段。但是这个方式用在大基因组上就非常没有效率了。我讲一下我的经历:

Step 1. 研究生阶段,为了实现给定一个基因组区间列表,快速从基因组里取序列,做了这么一件事情:把基因组按照染色体读到内存里,另外有一个辅助的Hash保存每个染色体的名字和长度。这个办法的麻烦之处在于每次都要把基因组load到内存里。其实后面还尝试了些优化,比如先统计Bed文件涉及到的chrom,只load这些。把ATCG读取成2进制,用2个bit表示。

Step 2. 后来接触到UCSC,开始研究Jimmy Kent 的UCSC格式,发现里面有个2bit格式正是用来解决这个问题的。其实我之前尝试的方案已经很接近2bit的方案了:把基因组转变成2bit,用00,01,10,11表示T,C,A和G。对于N的区间,额外用一个mask数据来记录。转换之后不是load在内存里,而是写在硬盘上。取序列时直接硬盘寻址。这个文件格式的好处是,存储文件经过2bit压缩,只有原来的1/4大小。这对于UCSC数据库是非常有意义的,毕竟在那时,硬盘还不是那么便宜。我在Python里写过一个类,专门读取2bit格式。这里把2bit文件的结构附上供有兴趣的人参考一下:

Step 3. 最后的解决方案是,也是现在最常用的解决方案,samtools。其实samtools做了一件非常简单的事情,简单到一说你就懂。通常基因组文件的存储形式为Fasta格式,就是一个标题行,跟着一长串的序列。为了显示方便,一般基因组文件的序列文件每行都是固定长,比如80bp(适应DOS屏幕显示)。既然有这个规律,我们就可以创建一个index文件,来记录这5个信息:(1)染色体名字,(2)染色体长度,(3)染色体的起始物理偏移地址,(4)每行序列的长度(eg. 60bp),(5)每行序列+换行符的长度 (Linux:61, Win:62)(windows为“\r\r",Linux为“\n",长度不同)。这样的话,给定任何一个染色体区间,我们就可以用这些信息计算出这个区间在硬盘上的存储位置,然后读取相应的序列。

对的,samtools 的faidx做的就是这个事情。hg19为例,hg19.fa.fai文件内容是这样的:
chr1    249250621       6       60      61
chr2    243199373       253404811       60      61
chr3    198022430       500657513       60      61
chr4    191154276       701980323       60      61
chr5    180915260       896320510       60      61
chr6    171115067       1080251031      60      61
chr7    159138663       1254218022      60      61


samtools的这个解决方案的好处是,基因组文件不需要进行格式转换。跟2bit相比,占用硬盘的空间大了些,但是考虑到现在硬盘不那么值钱了,加上samtools的普遍安装,这个解决方案是现在为止最popular的方案了。

有兴趣的同学可以写个简单的函数来实现下读取samtools Index好的基因组序列。Python代码大概是这样(未测试):
[Python] 纯文本查看 复制代码
import pandas
fai = pandas.read_table("hg19.fa.fai",index_col=0,header=None)
def getOffSet(chrom,start,end):    
    s,o,n,l = fai.ix[chrom,:]
    ostart = o + start/n*l+start%n
    oend   = o + end/n*l+end%n
    return ostart, oend

fh = open("hg19.fa")
for line in open('regions.bed'):
    chrom,start,end = line.rstrip().split('\t')
    ostart, oend = getOffSet(chrom,int(start),int(end))
    fh.seek(ostart,0) # move to start offset
    seq = "".join(fh.read(oend-ostart).split())
    print ">{0}:{1}-{2}\n{3}".format(chrom,start,end,seq)
fh.close()


继续向大家安利Python。一方面是我喜欢Python的编程思想和风格:一个问题有多个解决方案,我只选最好的那一个!另一方面,也是对我来讲更重要的方面,Python是一种只需要load几个有限的库,就可以做许多许多事情的语言,非常适合记忆力差,智商又不够用的我。Perl我也曾经学过两年,感觉就是学不会啊学不会,R也是在门外徘徊许久始终不得其门而入。对于BioPerl和BioConductor,我更是一进去就迷失,我看不懂啊看不懂。。。


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【课程预告】手把手教你入门生信——The Biostar Handbook
下一篇:R的bioconductor里面的stringDB包来做PPI分析
回复

使用道具 举报

0

主题

8

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2017-3-1 11:35:30 | 显示全部楼层
太棒了!终于直到.fai文件是干啥用的了!
我也很喜欢python,以后就可以尝试用这种方式来批量截取序列啦!
回复 支持 反对

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
 楼主| 发表于 2017-3-1 23:56:32 | 显示全部楼层
本帖最后由 tsznxx 于 2017-3-1 23:57 编辑
银色麦穗 发表于 2017-3-1 11:35
太棒了!终于直到.fai文件是干啥用的了!
我也很喜欢python,以后就可以尝试用这种方式来批量截取序列啦! ...

这个功能被整理到pysam里面了。Python讲究的是把任何好用的东西写成function,整理成class。然后在需要的时候只需要调用函数就好,像这样:
[Python] 纯文本查看 复制代码
pysam.faidx(infile) # create index
fh = pysam.Fastafile(infile) # open fasta file
for bed in BEDFILE:
    seq = fh.fetch(reference=chrom,start=start,end=stop) # fetch sequence
    # post analysis
fh.close()

而不是在shell里那样,把先文件做成bed文件,然后samtools读bed文件和Fasta文件,提取序列写到文件,最后在后续的script里读取序列文件,进行下一步处理。
pysam里面还整理很多很好的东西,有空了我会拿出来讲一下。
回复 支持 反对

使用道具 举报

0

主题

8

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2017-3-3 21:33:03 | 显示全部楼层
tsznxx 发表于 2017-3-1 23:56
这个功能被整理到pysam里面了。Python讲究的是把任何好用的东西写成function,整理成class。然后在需要的 ...

非常感谢!
回复 支持 反对

使用道具 举报

7

主题

15

帖子

116

积分

注册会员

Rank: 2

积分
116
QQ
发表于 2017-3-24 20:33:08 | 显示全部楼层
好像很厉害啊,学习一下。有问题再请教,嘿嘿
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 09:08 , Processed in 0.034937 second(s), 29 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.