了解到在seq中multi-hits的定义是 unique best hit, 即是read比对到了不同的位置，但是获得了同样的最高分数。而bowtie2默认的 >1 alignments 统计的是分数不低于 threshold 的 alignments 的数量。因此两种统计数字会出现偏差。
for FN in *.fa; do
echo $FN
GN=`echo $FN | cut -d'.' -f2`
bowtie2 -k 2 --end-to-end -f --large-index -x ~/*basename -U ${FN} -S `echo $FN | sed s/.fa//g`.bowtie2.${GN}.sam 2> `echo $FN | sed s/.fa//g`.bowtie2.${GN}.log &
done
(for i in *.bowtie2.*.sam; do
echo -e $i"\t"`gawk '!/^@/' $i | cut -f1,3,12 | gawk 'NR%2==1{ printf("%s\t%s\t%s\t", $1, $2, $3);} NR%2==0{ printf("%s\t%s\n", $2, $3);}' | sed 's/AS:i://g' | gawk '{T++; if($3==$5){C++;}} END{printf("%.2f%\n", C*100/T);}'`
done ) > multiple_hits.bowtie2
关于重定向
liunx 的file descriptor有三种 stdin、 stdout、 stderr 。其中重定位符>默认是 stdout，因此如果在指令出现错误的情况下，错误信息并不会被重定位至文件中。操作 stderr 的方法是 2> 。 如果想将 stdout、 stderr 重定位至同一个文件中，方法为command > file 2>&1 。2>&1在后面的原因是
command > file 2>&1
command > file
将标准输出重定向到file中， 2>&1 是标准错误拷贝了标准输出的行为，也就是同样被重定向到file中，最终结果就是标准输出和错误都被重定向到file中。
command 2>&1 >file
2>&1 标准错误拷贝了标准输出的行为，但此时标准输出还是在终端。>file 后输出才被重定向到file，但标准错误仍然保持在终端。
接上文
代码第一行的行为是使用bowtie进行比对，并将错误信息保存至日志文件中。第二行的行为是进入bowtie比对后的sam格式文件中将alignment的分数信息提取出来，两两比对是否相同。最后计算aligment综述与分数相同alignment数量的比值。这段代码的目的是某一次比对的 unique best hit 的比例。但是这样计算出来的结果可能并不正确。
在使用bowtie2进行比对的时候使用了参数-k 2。使用这个参数的原因是需要统计的量是multi-hits的数量，而只要只要找到一个以上，即2个，即可将这个read计入multi-hits中。
引用bowtie2 manual的说法
When -k is specified, however, bowtie2 behaves differently. Instead, it searches for at most <int> distinct, valid alignments for each read. The search terminates when it can't find more distinct valid alignments, or when it finds <int>, whichever happens first.
这说明即便使用了-k 2参数，也不能保证每个 read 都能有2个 alignment ，而这正是上面那段代码的基础。下面我进行了一些验证。
首先，我模拟了10000个随机的read。按照预想的情况-k 2参数应该保证在sam文件中，应该会有精确的20000行alignmen信息。而实际上alignment信息只有13528行，所以这个前提条件是不可能成立的。
另外一方面，两个alignment的所属的read名称应该相同。而实际上有很多单独的alignemnt。
那么使用这样可能会不成立的前提条件会导致什么结果呢。可以肯定的是，并不会有什么错误信息出现。对于模拟数据而言，由于没有加入点突变的错误率，所有read都是基因组的无错误复制，导致大量的alignment的分数都是 end-to-end 模式下的最高分0分，进而导致大量实际上并没有组成unique best hit的pair分数相同，进而导致被认定为multi-hits，进而导致比率虚高。
那么正确的比率应该如何计算呢。
- 确认两个alignment属于同一个read
- 确认两者的分数相等
修正后的比率有很大的降低。
关于bash脚本的执行过程
上面那段代码逻辑上可以分为两个部分，第一部分调用bowtie程序并将alignment结果保存在sam格式的文件中。第二部分代码读取sam文件中的内容，提取alignment信息统计multi-hits数量。执行顺序应为等待第一部分执行完毕后，等待sam文件生成完整，再执行第二部分代码。然而在实际执行过程时，几乎是瞬间提示不存在.sam文件，这说明并不是等待结束的顺序执行。查到一个命令 wait 在两部分代码之间加入这个命令可以过得理想顺序。
一个疑问。我之前写的一个脚本的执行顺序与此不同。有待检查。
网友评论