2

我正在尝试利用 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 中执行此操作，因为我担心循环内的循环会太慢。

4

1 回答 1

0

但是，这不会遍历所有文件，并且似乎卡在 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

注意：抛开[[条件中的变量，您应该注意引用可能包含一些空格的变量，例如文件名。否则它会破裂。

于 2020-12-31T13:26:28.317 回答