整数编码VCF文件的最优解

Question

我有解决问题的有效方法，但速度很慢。我很好奇推荐的加速方法，看看它能多快。这是一个示例输入文件

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   0|1:7,28:35:99:0|1:616504_T_C:787,0,177:616504   
Chr16 616504 T     C     BESC.1    0/0:48,0:48:99:.:.:0,114,1710:.                  
Chr16 616504 T     C     BESC.10   1|1:0,23:23:72:1|1:616504_T_C:1059,72,0:616504   
Chr16 616504 T     C     BESC.100  0/0:34,0:34:96:.:.:0,96,1440:.                   
Chr16 616504 T     C     BESC.1001 0/0:47,0:47:99:.:.:0,120,1800:.                  
Chr16 616504 T     C     BESC.1002 0/0:39,0:39:99:.:.:0,108,948:.

目标是从 value 列中取出第 1 个和第 3 个字符并将它们相加，然后输出一个类似的文件，其中值列替换为这个和。前两行的示例输出：

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   1   
Chr16 616504 T     C     BESC.1    0

这是我当前的解决方案，其中 STDIN 1 是输入文件名，STDIN 2 是输出文件名：

#!/bin/bash
i=0
len=$(cat  | wc -l)

touch 
while read -r line; do
    let "i++"
    geno=$(echo "$line" | awk '{ print  }'| cut -c1,3 | fold -w1 | paste -sd+ - | bc)

    echo "$line" | awk -v g="$geno" '{ =""; print [=12=] " " g}' >> 

    echo "Processed " "$i/$len"

done <

实际输入文件有1,707,993个条目。使用我的解决方案，这大约需要 4.5 小时来计算。理想情况下，我可以运行在一个小时内完成此操作，但我不确定这有多现实。谢谢！

Answer 1

这个循环会很慢，因为它在输入文件的每行中创建和销毁 6 个外部进程。

在 Bash 中，调用一个命令来处理整个文件通常比在该文件中每行调用一次命令要快得多。

例如，如果你想在 Awk 中进行文件处理，你可以这样做：

#!/bin/bash
set -eu

awk -v out_file="" 'NR == 1 {
    print > out_file;
}
NR != 1 {
     = substr(, 1, 1) +  substr(, 3, 1);
    print > out_file;
    print "Processed " NR;
}' ""

说明：

取第6个字段，取第一个和第三个字符求和。然后 re-assigns 他们回到字段 6。在 awk 中，这意味着 print 将发出原始输出，其中一个字段已更改。
打印重定向到 out_file。
发出进度指示器。

整数编码VCF文件的最优解

Optimal Solution for Integer Encoding VCF file

bash

awk

bioinformatics

vcf-variant-call-format