今天一直在尝试用python获取基因序列,原理嘛,就跟爬虫一样,但是巨慢,获取大量的基因序列的时候慎用!!,而且中途很容易出错。
然后在网上查找到,用bedtools真香
以下主要参考博客
bedtools批量提取基因组指定位置序列 - 简书 (jianshu.com)
我的是要从人类基因里面提取序列
首先我要获取到人类基因的fa文件
#下载人类基因组信息
wget url=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
#解压
tar zvfx chromFa.tar.gz
#将所有染色体信息做成一个fa文件
cat *.fa > hg19.fa
然后了解bed格式是什么样子滴~
BED文件必须的3列:
chrom - 染色体号; 例如,chr1,chrX。。。。。。。
chromStart - feature在染色体上起始位置. 从0开始算,染色体上第一个碱基位置标记为0。
chromEnd - feature在染色体上终止位置。染色体上前100个碱基片段的位置位置标记为:chromStart=0, chromEnd=100。 实际上,第100个碱基不属于当前片段中,当前片段的碱基应该是0-99。所以在BED文件中,起始位置从0开始,终止位置从1开始。
BED文件可选的9列:
4.name - BED行名,在基因组浏览器左边显示;
5.score - 在基因组浏览器中显示的灰度设定,值介于0-1000;
6.strand - 正负链标记. Either "." (=no strand) or "+" or "-".
7.thickStart - feature起始位置。绘制特征的起始位置(例如,基因显示中的起始密码子)。当没有这部分时,thickStart和thickEnd通常设置为chromStart位置。
8.thickEnd - feature编码终止位置
9.itemRgb - R,G,B (e.g. 255,0,0)值,当itemRgb 设置为 "On",BED的行会显示颜色
10.blockCount - blocks (exons)数目
11.blockSizes - blocks (exons)大小列表,逗号分隔,对应于blockCount
12.blockStarts -blocks (exons)起始位置列表,逗号分隔,对应于blockCount.;这个起始位置是与chromStart的一个相对位置。
我把文件处理成了这样,然后把后缀名改成了.bed(我只有染色体号,起始位点,终止位点)
最后,就可以两个文件都准备好了就可以开始bedtools了
$bedtools getfasta -fi 参考基因组.fa -bed 提取序列.bed -fo 输出文件.fa
除了转换成bed格式的文件实在自己电脑弄的,其他都在服务器上弄的