根据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 SeqIOdic={}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/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？相关推荐

gff文件_根据gff/gtf等注释文件取负链上的序列：先反向互补染色体再截取？还是先截取区间再反向互补序列？...
最近需要根据注释文件在基因组上截取序列,突然想到一个问题:对于下面这样在负链上的基因,我们是先将整条染色体反向互补再截取对应区间?还是先截取对应区间再反向互补呢? 首先亮出答案:先截取区间,再反向互补 ...
GTF基因注释文件详解
GFF和GTF是两种最常用的数据库注释格式,在信息分析中建库时除了需要fasta文件一般还会需要这两种文件,提取需要的信息进行注释. Cufflinks/Tophat 软件需要 GTF文件作为基因注释 ...
FX2N通讯监控记录， PLC程序/参数/注释/文件寄存器的写入和读出
//============================================================================================// // ...
gffread gffcompare 将gff与gtf格式的注释文件转换与合并
gffread gffcompare 将gff与gtf格式的注释文件转换与合并使用: (1)gffread 安装: conda install gffread -y 使用: mkdir gtf# 格 ...
「小技巧」如何让IGV更快的加载GTF和GFF注释文件
很简单,就下面3行命令 gff= (grep ^"#" $gff; grep -v ^"#" $gff | sort -k1,1 -k4,4n) | bgzip ...
根据gtf格式的基因注释文件得到人所有基因的染色体坐标
用bedtools对基因组片段区域进行基因注释根据gtf格式的基因注释文件得到人所有基因的染色体坐标选择的genecode内最早的Grch38版本(201408) v20是最早的hg38版本对应的 ...
【GO富集分析】GO注释文件爬取
GO数据库注释文件爬取爬取整体思路代码实现最近在做基因富集分析发现,很多非模式植物通过 clusterprofiler做富集分析都需要自备注释文件,这时我们需要GO的注释文件,需要自己整理,这里 ...
关于基因组注释文件GTF的解释
GTF文件的全称是gene transfer format,主要是对染色体上的基因进行标注.怎么理解呢,其实所谓的基因名,基因座等,都只是后来人们给一段DNA序列起的名字而已,还原到细胞中就是细胞核里 ...
构建index所需的参考基因组以及各种版本的注释文件
文章目录一.参考基因组 1. UCSC 2. ensemble 3. NCBI 4. gencode 二.基因组注释文件(GFF,GTF) 1. UCSC 2. ensemble 3. NCBI 4 ...

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

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

最新文章

热门文章