cd ~/DATA/Data_Tam_Methylation_19606WT_adeAB_adeIJ_craA

# ============================================================
# 0) 建库：核酸全库 + Acinetobacter 核酸子库 + 蛋白全库（blastx 用）
# ============================================================
makeblastdb -in ISfinder-sequences/IS.fna -dbtype nucl -title ISfinder 2>/dev/null || \
  echo "nucl DB 已存在，跳过"

awk '/^>/{keep = ($0 ~ /ISAba|Acinetobacter/)} keep' ISfinder-sequences/IS.fna > ISfinder-sequences/Acinetobacter_IS.fa
makeblastdb -in ISfinder-sequences/Acinetobacter_IS.fa -dbtype nucl -title Acinetobacter_IS 2>/dev/null || \
  echo "Acinetobacter nucl DB 已存在，跳过"

# 蛋白库（转座酶）——blastx 比 blastn 敏感得多，能找回 blastn 漏掉的 IS21/IS110
if [ -f ISfinder-sequences/IS.faa ]; then
  makeblastdb -in ISfinder-sequences/IS.faa -dbtype prot -title ISfinder_prot 2>/dev/null || \
    echo "prot DB 已存在，跳过"
else
  echo "⚠️ 未找到 IS.faa，将跳过 blastx 步骤"
fi

# ============================================================
# 1) 定义 4 个基因组（ISEScan 目录 → 对应 FASTA）
# ============================================================
REF_FA=$(find . -maxdepth 3 -name "CP059040.fasta" | head -1)

declare -A GENOMES=(
  [isescan_19606WT_polypolished]=unicycler-medaka_polished_genomes/19606WT_polypolished.fasta
  [isescan_19606_adeAB_polypolished]=unicycler-medaka_polished_genomes/19606_adeAB_polypolished.fasta
  [isescan_19606_adeIJ_polypolished]=unicycler-medaka_polished_genomes/19606_adeIJ_polypolished.fasta
  [isescan_CP059040]=$REF_FA
)

for fa in "${GENOMES[@]}"; do samtools faidx "$fa"; done

# ============================================================
# 2) 逐基因组：提取 IS → blastn → blastx → 合并命名
# ============================================================
for d in "${!GENOMES[@]}"; do
  fa=${GENOMES[$d]}
  gff=$(find $d -name "*.gff" ! -path "*proteome*")
  seqs=${d}_IS_seqs.fa
  > $seqs

  # -- 2a. 只提取带 family= 的行（即 8 个 IS 本体，排除 TIR）--
  grep -v "^#" "$gff" | awk -F'\t' '$9 ~ /family=/' | \
  while IFS=$'\t' read -r chr src typ start end score strand phase attr; do
    fam=$(echo "$attr" | grep -o 'family=[^;]*' | head -1 | cut -d= -f2)
    samtools faidx "$fa" "${chr}:${start}-${end}" | \
      awk -v h=">${chr}_${start}_${end}_${strand}_${fam}" 'NR==1{print h;next}{print}' >> $seqs
  done

  # -- 2b. blastn 全库 --
  blastn -query $seqs -db ISfinder-sequences/IS.fna \
    -evalue 1e-20 -max_target_seqs 1 -max_hsps 1 \
    -outfmt "6 qseqid sseqid pident length evalue bitscore stitle" \
    -out ${d}_IS_names_fullDB.tsv

  # -- 2c. blastn Acinetobacter 子库 --
  blastn -query $seqs -db ISfinder-sequences/Acinetobacter_IS.fa \
    -evalue 1e-20 -max_target_seqs 1 -max_hsps 1 \
    -outfmt "6 qseqid sseqid pident length evalue bitscore stitle" \
    -out ${d}_IS_names_Acinetobacter.tsv

  # -- 2d. blastx 蛋白库（更敏感）--
  if [ -f ISfinder-sequences/IS.faa ]; then
    blastx -query $seqs -db ISfinder-sequences/IS.faa \
      -evalue 1e-10 -max_target_seqs 3 -max_hsps 1 \
      -outfmt "6 qseqid sseqid pident length evalue bitscore stitle" \
      -out ${d}_IS_names_blastx.tsv
  fi

  # -- 2e. 三级合并命名：blastn全库 → blastx蛋白库 → 标记无近缘 --
  # 第 1 级：blastn 全库命中（每个 query 取 1 行）
  awk -F'\t' -v OFS='\t' '!seen[$1]++ {print $0, "blastn_fullDB"}' \
      ${d}_IS_names_fullDB.tsv > ${d}_IS_final_names.tsv

  # 第 2 级：blastn 未命名的 query，用 blastx 最佳命中补充
  if [ -f ${d}_IS_names_blastx.tsv ]; then
    awk -F'\t' -v OFS='\t' 'NR==FNR{named[$1]=1; next} !($1 in named) && !seen[$1]++ {print $0, "blastx_protDB"}' \
        ${d}_IS_names_fullDB.tsv ${d}_IS_names_blastx.tsv >> ${d}_IS_final_names.tsv
  fi

  # 第 3 级：仍无命中的 IS，标记为 no_close_relative（保证 8 行齐全）
  grep "^>" $seqs | sed 's/^>//' | while read q; do
    if ! cut -f1 ${d}_IS_final_names.tsv | grep -qxF "$q"; then
      printf '%s\t-\t-\t-\t-\t-\t-\t-\tno_close_relative_in_ISfinder\n' "$q" >> ${d}_IS_final_names.tsv
    fi
  done

  echo "==================== $d ===================="
  # qseqid=坐标_链_ISEScan家族 | sseqid/stitle=ISfinder命名 | 末列=命名方法
  column -t -s $'\t' ${d}_IS_final_names.tsv | cut -c1-180
done