用bedtools從基因組序列提取gene部分

參考生信寶典教程: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
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容