通常我们碰到的是，根据gff文件提取cDNA序列，这个比较容易满足，但是有时候我们需要构建native promoter 载体的时候，就需要提取cDNA前面的序列，真核生物大约是3k，原核大约是0.5k。这个时候就需要稍微绕一下弯，但是思路还是基本一致的。
1.准备的东西
工具：samtools,bedtools,bedops，都可以通过conda一键安装。
序列文件：基因组文件fasta，注释文件gff。
2.提取
2.1 提取gff文件的所有基因位置,并转换成bed格式
awk '{if($3~/^gene$/)print}' file.gff > genes.gff
gff2bed <genes.gff> genes.bed
2.2 计算染色体长度
samtools faidx file.fasta
cut -f 1,2 file.fa.fai > sizes.chr
2.3 创建包含promoter位置的bed文件
bedtools flank -i genes.bed -g sizes.chr -l 3000 -r 0 -s > promoters.bed
# -l 基因起始位置前多少bp
# -r 基因后多少bp
# -s
#注意，这里生成的bed文件缺少基因id，手动添加即可，记住添加到第四列
2.4 根据promoter的位置信息，在基因组序列中抓取promoter的序列
bedtools getfasta -s -fi file.fa -bed promoters.bed -fo promoters.fa -name
网友评论