从单个字符串有效地计算统计数据（bowtie2 的 bam 文件）

Question

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

我们从这样的字符串开始

    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  | 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/ )
    fi
    }

然而，这个解决方案是错误的，因为它将 ^TCTAAG 等删除视为不匹配。

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

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

感谢您的考虑

Answer 1

使用二进制数据（bam 文件）输入和输出要快得多，但是你有一个 sam 文件（平面版本）
正如@cel 之前已经评论过的，用编译语言（C、C++ 等）编写的程序总是比脚本（ej、bash、awk 等）快

我在awk中展示了一个脚本，供您尝试

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

你得到

96003934 98.31

说明

第 18 列：MD:Z:54G53T9^TACA11

匹配 = 54+53+9 = 116

不匹配 = count_letter(54G53T9) = 2

id = 116 / (116+2) = 0.9830508474576272

Answer 2

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

awk '{
 match([=10=], /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", , cmatch*100/(cmatch+nmatch-1));
}'

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

从单个字符串有效地计算统计数据（bowtie2 的 bam 文件）

Efficiently calculate statistics from single string (bam file by bowtie2)

optimization

parsing

bioinformatics