統計測序數據reads數和堿基數的幾種方法

2024年2月6日 21点热度 0人点赞

很簡單的問題,卻被常常問起。記錄一個帖子。文末有福利。

手動寫一個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

我們的論壇建好了,測試了一段時間,用起來也正常了。在上次的調查中,大部分老師也提議通過論壇交流,一來問答清晰,二來便於追蹤記錄,後續遇到問題直接搜索下就可以有答案,網址是http://www.ehbio.com/Esx/ (微信掃碼就可以登錄)