fastx_toolkit由一系列的命令组成,每个命令提供一个实用的小功能。在使用时需要注意以下几点
- 不支持压缩格式的输入文件
- 不允许序列中存在N碱基,这样的序列会自动去除
- 可视化命令依赖gunplot软件和perl的GD模块
- 默认情况下认为fastq文件的碱基编码格式为phred64,对于phred33编码的fastq文件,需要添加参数-Q 33
可以利用以下代码判断fastq文件编码是phred33还是64
fastq格式文件及phred33的判断
#!/bin/bash
# phed.sh
# Usage: ./phed.sh <fastq_file>
# 确保有一个文件名作为输入
if [ -z "$1" ]
then
echo "Usage: $0 <fastq_file>"
exit 1
fi
# 分析前1000行中的质量得分
zless "$1" | head -n 1000 | awk '{if(NR%4==0) printf("%s",$0);}' \
| od -A n -t u1 -v \
| awk 'BEGIN{min=100;max=0;} \
{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END \
{if(max<=126 && min<59) print "Phred33"; \
else if(max>73 && min>=64) print "Phred64"; \
else if(min>=59 && min<64 && max>73) print "Solexa64"; \
else print "Unknown score encoding"}'