统计测序数据reads数和碱基数的几种方法

很简单的问题，却被常常问起。记录一个帖子。文末有福利。

手动写一个FASTQ格式的测试数据

cat <<END >sample.fq
@ESX1
CAGGAGGAGTACGTGTTTTTTTTTTGCAGTACTGTACGGCGCAGTAC
+
FFFFFFFFFFFFFFEEFFFFFFFFFFFFFFFFFFFFFEEEFFFFFFF
@ESX2
CAGGAGGAGTACGTGTTTTATTTTTGCAGTACTGTACGGCGCAGTAC
+
FFFFFFFFFFFFFFEEFFFFFFFFFFFFFFFFFFFFFEEEFFFFFFF
@ESX3
CAGGAGGAGTACGTGTTTTTTTTTTGCAGTACTGTACGGCGCAGTAC
+
FFFFFFFFFFFFFFEEFFFFFFFFFFFFFFFFFFFFFEEEFFFFFFF
END

利用seqkit统计

更详细的介绍和安装见推文seqkit：序列梳理神器-统计、格式转换、长度筛选、质量值转换、翻译、反向互补、抽样、去重、滑窗、拆分等30项全能。

可以同时统计单个或多个fastq文件，结果输出为表格形式

seqkit stat sample.fq
# 结果如下
# num_seq：总序列数
# sum_len: 总碱基数
file       format  type  num_seqs  sum_len  min_len  avg_len  max_len
sample.fq  FASTQ   DNA          3      141       47       47       47

# 统计多个文件
seqkit stat sample.fq sample.fq
file       format  type  num_seqs  sum_len  min_len  avg_len  max_len
sample.fq  FASTQ   DNA          3      141       47       47       47
sample.fq  FASTQ   DNA          3      141       47       47       47

# 统计多个压缩文件
seqkit stat *.fq.gz

用Linux命令统计

awk的介绍见常用和不太常用的awk命令

# 统计单个文件
# awk运算
# %取余数
# 为什么除以4，又除以1000000？cat sample.fq | awk 'BEGIN{OFS="\t"}{if(FNR%4==0) base+=length}END{print FNR/4/1000000 " million", base/10^9 "G";}'
# 3e-06 million 1.41e-07 G

# 统计多个文件
for i in *.fq; do 
  cat sample.fq | awk -v name=${i} 'BEGIN{OFS="\t"}{if(FNR%4==0) base+=length}END{print name, FNR/4/1000000 " million", base/10^9 " G";}'
done

# sample.fq       3e-06 million   1.41e-07 G

# 统计多个压缩文件
for i in *.fq.gz; do 
  zcat sample.fq.gz | awk -v name=${i} 'BEGIN{OFS="\t"}{if(FNR%4==0) base+=length}END{print name, FNR/4/1000000 " million", base/10^9 " G";}'
done