整数编码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列中取出第一个和第三个字符并对它们求和，然后输出一个类似的文件，其中值列替换为该总和。前两行的示例输出：

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 $1 | wc -l)

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

    echo "$line" | awk -v g="$geno" '{ $6=""; print $0 " " g}' >> $2

    echo "Processed " "$i/$len"

done < $1

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

Answer 1

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

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

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

#!/bin/bash
set -eu

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

解释：

1. 它采用第 6 个字段，并采用第一个和第三个字符并对它们求和。然后它将它们重新分配回字段 6。在 awk 中，这意味着 aprint将发出更改了第一个字段的原始输出。
2. 打印被重定向到out_file.
3. 发出进度指示器。