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

2021-10-28  本文已影响0人  学生信的大叔

参考生信宝典教程: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
上一篇 下一篇

猜你喜欢

热点阅读