参考生信宝典教程：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
网友评论