用bedtools从基因组序列提取gene部分

参考生信宝典教程:https://blog.csdn.net/qazplm12_3/article/details/118316221

提取基因序列的操作也类似于提取启动子序列。这里要注意GFF文件的序列位置是从1开始,而bed文件的位置是从0开始,前闭后开,所以要对序列的起始位置进行-1的操作。
修改了原代码提取的列部分,用 基因ID列 $10 作为提取基因序列后的序列ID。

type="gene"
sed 's/"/\t/g' GCA_003024255.1_Red5_PS1_1.69.0_genomic.gtf |awk -v type="${type}" 'BEGIN{OFS=FS="\t"}{if($3==type) {print $1,$4-1,$5,$10,".",$7}}' > red5.geneid.bed
bedtools getfasta -name  -fi GCA_003024255.1_Red5_PS1_1.69.0_genomic.fna -bed red5.geneid.bed >Red5.geneid.fa
#核对基因数目
#cat GCA_003024255.1_Red5_PS1_1.69.0_genomic.gtf |awk '{if(($3 == "gene")) print  }' | wc -l
#cat Red5.geneid.fa |grep ">" |wc -l
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容