在循环中循环 vcftools bash

问题描述

我正在尝试利用 vcftools 包来计算堰和科克汉姆的 fst。我想在第一个实例中遍历两对种群，然后在 1000 Genomes 项目的所有变体中循环这些种群：每个染色体包含一个单独的 vcf 文件。例如，对于 pop1 与 pop2，对于 pop3 与 pop4，计算染色体 1-10 的 fst。每个种群文件，例如 LWKfile 包含属于该种群的个体列表。

我已经尝试过：

for population in LWK_GBR YRI_FIN; do

firstpop=$(echo $population | cut -d '_' -f1)
secondpop=$(echo $population | cut -d '_' -f2)

for filename in *.vcf.gz; do

vcftools --gzvcf ${filename} \
--weir-fst-pop /outdir/${firstpop}file \
--weir-fst-pop /outdir/${secondpop}file \
--out /out/${population}_${filename}

done

done  

然而，这不会遍历所有文件，似乎卡在 10 号染色体上。有没有更有效的方法在 bash 中执行此操作，因为我担心循环内的循环会太慢。

解决方法

然而，这并没有遍历所有文件，似乎卡在了 10 号染色体上。

我担心循环内的循环会太慢。

您确定是 for filename in *.vcf.gz 太慢而无法循环所有文件吗？

尝试在 echo 之前放一个 vcftools 以查看它是否仍然卡住。

您需要确定什么需要花费太多时间才能做出正确的选择。

例如如果它是 vcftools 可能你不需要等待这个命令的结束并考虑做一些异步处理。

如果一个循环的文件太多，你也应该考虑做一些并行处理。

此外，您似乎对所有 .vcf.gz 文件重复了两次循环。反转你的两个循环可能会快一点。

以下是使用 bash 进行并行和异步处理的示例：

#!/bin/bash

MAX_PARALLEL_PIDS=4 # adjust regarding your own machin capacity (cpu available,etc... it could be dynamically calculated)

declare -a POPS
declare -a PIDS

POPS=("LWK_GBR" "YRI_FIN")

# your heavy treatment in a function
process() {
  pop="${1}"
  filename="${2}"
  firstpop="${pop%%_*}" # no need to call an external program here
  secondpop="${pop#*_}" # same here

  vcftools --gzvcf "${filename}" \
     --weir-fst-pop "/outdir/${firstpop}file" \
     --weir-fst-pop "/outdir/${secondpop}file" \
     --out "/out/${pop}_${filename}"
}

# a function which is usefull to wait all process when your "thread pool" reached its limits
wait_for_pids() {
  for pid in "${PIDS[@]}"; do
    [[ $pid =~ ^[0-9]+ ]] && wait $pid
  done
  unset PIDS
}

i=0
for filename in *.vcf.gz; do
 if [[ $i -ge $MAX_PARALLEL_PIDS ]]; then
   i=0
   wait_for_pids
 fi

 for population in "${POPS[@]}"; do
   process "${population}" "${filename}" & # You won't wait for the end here
   PIDS[$i]=$!
   (( i++ ))
 done
done

# at the end wait for the remaining processes
wait_for_pids

N.B: 撇开 [[ 条件中的变量不谈，你应该注意引用可以包含一些空格的变量，尤其是文件名等。否则会坏掉。