写在前面
通常我们碰到的是,根据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
2.4 根据promoter的位置信息,在基因组序列中抓取promoter的序列
bedtools getfasta -s -fi file.fa -bed promoters.bed -fo promoters.fa -name