gff文件_根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？...

最近需要根据注释文件在基因组上截取序列，突然想到一个问题：对于下面这样在负链上的基因，我们是先将整条染色体反向互补再截取对应区间？还是先截取对应区间再反向互补呢？

首先亮出答案：先截取区间，再反向互补。比如上面的ALK基因，先截取chr2上的29415640-30144452区间，再反向互补即得到ALK基因的碱基序列。

验证过程：

方法一：bedtools工具包可以根据bed文件提取区间序列，这里以上面的ALK基因为例试一下：

1.首先创建一个bed文件命名为AKL.bed，其中文件第6列可以指定正负链

2.然后运行bedtools，得到ALK.fa序列文件，注意要添加-s选项

bedtools getfasta -fi ucsc.hg19.fasta -bed ALK.bed -s -fo ALK.fa

3.比较NCBI上给出的ALK序列和bedtools的运行结果，发现两者一致，说明bedtools工具截取负链上的基因准确，如果不想深究是哪种截取方法，后面直接用这个工具即可

方法二：将两种不同的截取方法截取的序列和NCBI上给出的ALK序列比较，即可知道哪种截取方式是正确的：

1.利用python编写一个用于获取反向互补序列的子函数

def rev(seq):        base_trans = {'A':'T', 'C':'G', 'T':'A', 'G':'C', 'N':'N', 'a':'t', 'c':'g', 't':'a', 'g':'c', 'n':'n'}        rev_seq = list(reversed(seq))        rev_seq_list = [base_trans[k] for k in rev_seq]        rev_seq = ''.join(rev_seq_list)        return(rev_seq)

2.利用python编写一个获取基因组所有序列的子函数

def cget_fadic(faPath):        from Bio import SeqIO        dic={}        for i in SeqIO.parse(open(faPath),'fasta'):                dic[str(i.id)]=str(i.seq)        return dic

3.截取序列并和NCBI比较

mydic=cget_fadic("ucsc.hg19.fasta")chr2_seq=mydic['chr2']seq1=rev(chr2_seq[29415640:30144452])seq2=rev(chr2_seq)[29415640:30144452]out1=open("ALK_seq1.fa",'w')out2=open("ALK_seq2.fa",'w')out1.writelines(">rev(chr2:29415640:30144452)\n"+seq1)out2.writelines(">rev(chr2):29415640:30144452\n"+seq2)out1.close()out2.close()

结果不难看出，seq1的序列和ALK序列一致，正确的序列是通过先截取再取反向互补的方式得到的。

gff文件_根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？...相关推荐

根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？
最近需要根据注释文件在基因组上截取序列,突然想到一个问题:对于下面这样在负链上的基因,我们是先将整条染色体反向互补再截取对应区间?还是先截取对应区间再反向互补呢? 首先亮出答案:先截取区间,再反向互补 ...
gff文件_[Py005] gff文件处理1
根据第3列的type,提取mRNA及相应exon的信息思路: 每次读取一行,提取到mRNA特征值后,写出该行: 判断下一行是否具有mRNA或exon特征值,如果有的话,递归自动判断下下一行 ...
GTF基因注释文件详解
GFF和GTF是两种最常用的数据库注释格式,在信息分析中建库时除了需要fasta文件一般还会需要这两种文件,提取需要的信息进行注释. Cufflinks/Tophat 软件需要 GTF文件作为基因注释 ...
用cmd运行python文件_怎么用cmd运行python文件
Layout Go工程项目的整体组织首先我们看一下整个 Go 工程是怎么组织起来的. 很多同事都在用 GitLab 的,GitLab 的一个 group 里面可以创建很多 project.如果我们进 ...
python运行系统找不到指定文件_系统找不到指定文件_系统找不到指定的文件_python 系统找不到指定的文件 - 云+社区 - 腾讯云...
广告关闭腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 我正在构建一个依赖于另一个库的特定库,当我这样做时,我收到以下警告消息:&quo ...
vbscript 写入文件_使用VBScript静默删除空文件夹
vbscript 写入文件 I decided to write this article after I saw an interesting question here on Experts Ex ...
计算机查找文件的速度,如何快速搜索文件_怎么加快电脑里的文件搜索速度
2019-11-01 10:40:56 浏览量:261 现在电脑的硬盘容量越来越大了,里面的内容也随之变多,所以当我们想要找到某个文件的时候也会变得很麻烦.特别是有些时间比较久了的,更加困难.那有网友 ...
python怎么打开ipynb文件_如何优雅的打开.ipynb文件
目前在windows打开.ipynb的文件的方法,网上几乎就只有一种,在cmd下 > jupyter notebook 这个方法只是方便你新建notebook的时候.如果你想再打开它,当你优雅的 ...
python多线程读取文件夹下的文件_是否可以使用python多线程从文件夹数读取文件数，并处理这些文件以获得组合结果？...
我认为学习使用线程的最简单方法是在concurrent.futures模块中使用ThreadPoolExecutor类,因为它比通常的同步for循环多了几行.尤其是在Python3中,但这可以适用于P ...

gff文件_根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？...

gff文件_根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？...相关推荐

最新文章

热门文章