#可以将测序基因组与已知基因集直接进行 blast 比对，然后将比对上的序列提取出来作为基因。
#基因组： NC_009648.1 https://www.ncbi.nlm.nih.gov/nuccore/NC_009648.1/ 
#质粒： NC_009649.1 https://www.ncbi.nlm.nih.gov/nuccore/NC_009649 
efetch -db nuccore -format fasta -id NC_009648.1 > MGH78578.fasta 
efetch -db nuccore -format fasta -id NC_009649 >>MGH78578.fasta

#下载参考序列基因集
#https://www.ncbi.nlm.nih.gov/genome/?term=NC_009648
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/240/185/GCF_000240185.1_ASM24018v2/GCF_000240185.1_ASM24018v2_protein.faa.gz
#建立索引
makeblastdb -in GCF_000240185.1_ASM24018v2_protein.faa -dbtype prot -parse_seqids -out GCF_000240185.1_ASM24018v2_protein.faa
#blastx 比对 blastx 
#核酸序列对蛋白库的比对，先将核酸序列翻译成蛋白序列（根据相位可以翻译为6种可能的蛋白序列），然后再与蛋白库做比对。
blastx -query MGH78578.fasta -out blast.out -db GCF_000240185.1_ASM24018v2_protein.faa -outfmt 6 -evalue 1e-5
#提取比对区域，生成 bed 文件
awk '{if ($7 < $8) print $1"\t"$7-1"\t"$8;else print $1"\t"$8-1"\t"$7}' blast.out >gene.bed
#根据比对位点，提取序列
seqkit subseq --bed gene.bed MGH78578.fasta >MGH78578_gene.ffn