parsing - 从单个字符串有效地计算统计数据(bowtie2 的 bam 文件)

Question

我的目标是有效地将包含短 DNA 测序读取的 bowtie2 映射的 bam 文件转换为包含映射开始和同一性百分比的简单表格。我即将完成此任务，但是我的解决方案非常慢并且无法处理重要的异常。我将逐步说明情况:

我们从这样的字符串开始

    FCC5G2YACXX:5:1101:1224:2059#NNNNNNNN   97  genome  96003934    24  118M4D11M   =   96004135    0   GCA....ACG  P\..GW^EO   AS:i:-28    XN:i:0  XM:i:2  XO:i:1  XG:i:4  NM:i:6  MD:Z:54G53T9^TACA11 YT:Z:UP

我们只取第四列和 MD:Z:54G53T9^TACA11 部分，后者表示匹配(数字)和不匹配(字母)。除非字母前面有“^”，否则字母是引用序列中的缺失。

然后我将其计算成这样的字符串 96003934 98.00

其中第一列与原始输入的第四列相同，第二列是匹配项，除以匹配项和不匹配项的总和。

因为我更喜欢 bash/zsh 脚本，所以我做了以下操作

    if_sam2tab() {
    tab=$(echo $1 | rev | cut -d '.' -f 2- | rev)
    if      [ ! -e ./bowtie_results/$tab.tab ]
    then    echo -e "rstart\tmatch\tmismatch" > ./bowtie_results/$tab.tab
            while   read -r l
            do      rstart=$(echo $l | cut -f 4 -d " "      )
                    t=$(echo $l | grep -o 'MD:Z:[0-9A-Z]*' )
                    match=$(echo ${t: 5} | tr '[a-zA-Z]' '+' | bc )
                    mismatch=$(echo ${t: 5} | tr -d '[0-9]\n' | wc -c )
                    sum=$(echo "$match + $mismatch" | bc )
                    id=$(echo "scale=2; $match / $sum *100"         | bc  )
                    echo -e "$rstart\t$id" >> ./bowtie_results/$tab.tab
            done < <(grep 'MD:Z' ./bowtie_results/$1 )
    fi
    }

然而，这个解决方案是错误的，因为它认为诸如 ^TCTAAG 之类的删除是不匹配的。

其次，它对我来说似乎非常慢，即使我在六个 cpu 上并行运行其中五个函数也是如此。

概括地说，我的目标是有效地将带有映射信息的字符串转换为标识百分比。

感谢您的关注

Answer 1

Jose 对先前回答的补充:MD 字符串不一定是第 18 列。因此，我在 Jose 的 awk 脚本中添加了一行来查找 MD 字符串，而不管行中的位置。毫无疑问，在这个过程中添加一个正则表达式步骤会减慢整个过程。

awk '{
 match($0, /MD:Z:[0-9A-Z\^]*/,m );
 split(m[0],v,/[\^:]/);
 nmatch = split(v[3],vmatch, /[^0-9]/);
 cmatch=0;
 for(i=1; i<=nmatch; i++) cmatch+=vmatch[i];
 printf("%s"OFS"%.2f\n", $4, cmatch*100/(cmatch+nmatch-1));
}' 

我选择在这个 awk 脚本中直接通过管道传输 bowtie 输出，而不将其保存在磁盘上。在我的有限测试中，我没有发现磁盘 i/o 受到限制，也没有发现将 bowtie2 输出直接输送到 awk 会显着降低同时运行时的整体性能。