如何不使用epost，用Entrez Direct批量查询核苷酸版本号并解决循环问题

我从NCBI Blast下载了核苷酸比对结果表（采用核苷酸集合数据库+megablast程序），先用awk按版本登录号排序：

awk -F "\t" 'NF>1{print}' unsorted_input.txt | sort -k2 > sorted_output.txt

接着用Entrez Direct通过版本登录号提取每条比对的目标物种：

awk -F "\t" 'NF>1{print $2}' unsorted_input.txt | epost -db nucleotide | efetch -format docsum | xtract -pattern DocumentSummary -element Organism | sort | paste sorted_output.txt - > final_output.txt

但这个命令只能提取部分物种数据，那些epost处理不了的条目，单独用esearch却能查询成功：

esearch -db nucleotide -query "accession_version_identifier" | efetch -format docsum | xtract -pattern DocumentSummary -element Organism

于是尝试用循环批量处理，提取每行第二列的版本登录号对应的物种：

while IFS=$'\t' read -r -a myArray
do
 esearch -db nucleotide -query "${myArray[1]}" | efetch -format docsum | xtract -pattern DocumentSummary -element Organism > "output.txt"
done < input.txt

但这个循环只返回第一行的物种。现在有两个问题：

1. 怎么修改循环，让所有行的物种都保存到同一个文件里？
2. 能不能不用epost，直接用Entrez Direct批量查询多个核苷酸版本登录号？

输入文件前几行（制表符分隔）：

ce1e013e-c4c5-47f9-b041-521ee293c4f0    AB002282.1  91.217  649 24  22  41  676 8   636 0.0 854
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  84.615  117 9   6   118 228 17668   17781   5.16e-19    108
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  84.615  117 9   6   118 228 20740   20853   5.16e-19    108
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  84.615  117 9   6   118 228 23812   23925   5.16e-19    108
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  84.615  117 9   6   118 228 26884   26997   5.16e-19    108
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  84.615  117 9   6   118 228 29956   30069   5.16e-19    108
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  85.345  116 9   6   118 228 33027   33139   1.11e-20    113
c10d7882-cc00-4ee2-8643-9b27fef66e83    AB828191.1  87.000  100 7   5   132 228 14613   14709   5.16e-19    108
8e8ac3f3-63f6-4519-ad25-287a25169f87    AB850654.1  88.262  4660    175 260 16  4401    103840  108401  0.0 5232
c4233926-9f23-46c4-bc4d-5702f47885bd    AB850654.1  89.958  4272    119 235 1   4042    104203  108394  0.0 5227
876d8f20-9d36-4207-8754-0924d99a6c46    AC019188.6  91.855  221 4   7   3   210 78509   78290   1.39e-75    296

1. 修改循环实现全量提取并保存到同一文件 #

原来的循环用>重定向会每次覆盖output.txt，改成>>追加模式即可。另外建议加一些防御性处理，避免空行或异常请求触发错误，同时加延迟避免NCBI的请求频率限制：

# 先清空输出文件，避免残留旧数据
> output.txt
# 循环处理每行，兼容文件末尾无换行的情况
while IFS=$'\t' read -r -a myArray || [ -n "${myArray[0]}" ]
do
  # 跳过空行或第二列无内容的行
  if [ -n "${myArray[1]}" ]; then
    esearch -db nucleotide -query "${myArray[1]}" | efetch -format docsum | xtract -pattern DocumentSummary -element Organism >> output.txt
    # 延迟1秒，降低触发NCBI限流的概率
    sleep 1
  fi
done < input.txt

如果需要把物种信息和原行内容对应起来，直接在循环里拼接后写入：

> final_output.txt
while IFS=$'\t' read -r line || [ -n "$line" ]
do
  acc=$(echo "$line" | awk -F "\t" '{print $2}')
  if [ -n "$acc" ]; then
    organism=$(esearch -db nucleotide -query "$acc" | efetch -format docsum | xtract -pattern DocumentSummary -element Organism)
    echo -e "$line\t$organism" >> final_output.txt
    sleep 1
  fi
done < input.txt

2. 不用epost，直接批量查询多个登录号 #

可以用esearch直接批量提交登录号，把多个登录号用逗号分隔即可。如果登录号数量多，建议分批次处理，避免触发NCBI的请求限制：

方法一：小批量一次性查询 #

先提取去重后的登录号，然后一次性查询，再和原表关联：

# 提取去重的登录号
awk -F "\t" 'NF>1{print $2}' input.txt | sort -u > unique_acc.txt
# 把登录号拼成逗号分隔的字符串
acc_list=$(paste -s -d ',' unique_acc.txt)
# 批量查询并提取登录号+物种对应关系
esearch -db nucleotide -query "$acc_list" | efetch -format docsum | xtract -pattern DocumentSummary -element AccessionVersion Organism > acc_organism.txt
# 关联原表和物种表（需要先对两个表按登录号排序）
sort -k2 input.txt > sorted_input.txt
sort -k1 acc_organism.txt > sorted_acc_organism.txt
join -1 2 -2 1 sorted_input.txt sorted_acc_organism.txt > final_output.txt

方法二：大数量分批次查询 #

如果登录号数量超过NCBI单次请求限制（一般建议不超过200个），可以拆分后逐个处理：

awk -F "\t" 'NF>1{print $2}' input.txt | sort -u > unique_acc.txt
# 每100个登录号分一个文件
split -l 100 unique_acc.txt acc_batch_
# 循环处理每个批次
> all_acc_organism.txt
for batch in acc_batch_*
do
  acc_list=$(paste -s -d ',' "$batch")
  esearch -db nucleotide -query "$acc_list" | efetch -format docsum | xtract -pattern DocumentSummary -element AccessionVersion Organism >> all_acc_organism.txt
  # 批次间延迟2秒，进一步降低限流风险
  sleep 2
done
# 关联原表和物种表
sort -k2 input.txt > sorted_input.txt
sort -k1 all_acc_organism.txt > sorted_all_acc_organism.txt
join -1 2 -2 1 sorted_input.txt sorted_all_acc_organism.txt > final_output.txt

内容的提问来源于stack exchange，提问作者Loggy01