參考生信寶典教程:https://blog.csdn.net/qazplm12_3/article/details/118316221
提取基因序列的操作也類似于提取啟動子序列。這里要注意GFF文件的序列位置是從1開始,而bed文件的位置是從0開始,前閉后開,所以要對序列的起始位置進(jìn)行-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
#核對基因數(shù)目
#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