一、fastq 文件格式
代码语言:javascript复制@DJB775P1:248:D0MDGACXX:7:1202:12362:49613
TGCTTACTCTGCGTTGATACCACTGCTTAGATCGGAAGAGCACACGTCTGAA
JJJJJIIJJJJJJHIHHHGHFFFFFFCEEEEEDBD?DDDDDDBDDDABDDCA
@DJB775P1:248:D0MDGACXX:7:1202:12782:49716
CTCTGCGTTGATACCACTGCTTACTCTGCGTTGATACCACTGCTTAGATCGG
IIIIIIIIIIIIIIIHHHHHHFFFFFFEECCCCBCECCCCCCCCCCCCCCCC
第一行:以‘@’开头,是这一条 read 的名字,这个字符串是根据测序时的状态信息转换过来的,中间不会有空格,它是每一条 read 的唯一标识符,同一份 FASTQ 文件中不会重复出现,甚至不同的 FASTQ 文件里也不会有重复;
第二行:测序 read 的序列,由 A,C,G,T 和 N 这五种字母构成,这也是我们真正关心的DNA 序列,N 代表的是测序时那些无法被识别出来的碱基;
第三行:以‘ ’开头,在旧版的 FASTQ 文件中会直接重复第一行的信息,但现在一般什么也不加(节省存储空间);
第四行:测序 read 的质量值,这个和第二行的碱基信息一样重要,它描述的是每个测序碱基的可靠程度,用 ASCII 码表示。
二、质量值体系
phred 质量值体系
illumina 测序质量值体系
从表中可以看到下限有 33 和 64 两个值,我们把加 33 的的质量值体系称之为 Phred33,加64 的称之为 Phred64(Solexa 的除外,它叫 Selexa64)。不过,现在一般都是使用 Phred33这个体系,而且 33 也恰好是 ASCII 的第一个可见字符('!')
文件格式:https://genome.ucsc.edu/FAQ/FAQformat.html#format1
三、fastq 格式文件处理
代码语言:javascript复制1 压缩与解压缩
解压缩
gunzip illumina_1.fastq.gz
gzip -d illumina_2.fastq.gz
压缩
gzip illumina_1.fastq
gzip illumina_2.fastq
2 fastq 文件统计
seqkit stats illumina_1.fastq.gz illumina_2.fastq.gz
3 统计 fastq 文件每条序列 ATCG 四种碱基组成以及质量值分布
seqtk comp illumina_1.fastq.gz illumina_2.fastq.gz
4 ATCG 以及质量值分布
seqtk fqchk illumina_1.fastq.gz
seqtk fqchk illumina_2.fastq.gz
57 交叉合并 pairend 文件
seqtk mergepe illumina_1.fastq.gz illumina_2.fastq.gz >merge.fastq
6 过滤短的序列
过滤小于 150bp 序列,并压缩输出
seqkit seq -m 150 nanopore.fastq.gz | gzip - >filter_150.fq.gz
保留小于 150bp 序列
seqkit seq -M 150 nanopore.fastq.gz
7 转换为列表格式
seqkit fx2tab nanopore.fastq.gz
8 分别统计每一条序列长度
seqkit fx2tab nanopore.fastq.gz |awk -F "t" '{print length($2) }'
9 质量值转换
将 illumina 1.8 转换为 1.5
seqkit convert --to Illumina-1.5 illumina_1.fastq.gz |head -4
将 illumina 1.5 转换为 1.8,什么都不加就是转换为 1.8
seqkit convert illumina_illmina1.5.gz
10 排序,按照长度
seqkit sort -l -r nanopore.fastq.gz
11 #seqkit 抽样,按照百分比
seqkit sample -p 0.1 illumina_1.fastq.gz
12 seqkit 抽样,按照条数
seqkit sample -n 1000 illumina_1.fastq.gz
13 拆分数据
seqkit split2 -1 illumina_1.fastq.gz -2 illumina_2.fastq.gz -p 2 -f
14 转换为 fasta
seqkit 工具
seqkit fq2fa nanopore.fastq.gz >nanopore.fasta
15 只输出 20 行 ID
seqkit seq -n -i nanopore.fastq.gz |head -20 >id.list
16 提取序列
seqkit grep -f id.list nanopore.fastq.gz
17 截取头尾
seqtk trimfq -b 15 -e 15 -Q illumina_1.fastq.gz
17 修改 reads ID
seqkit replace -p "SRR8494939.sra" -r 'reads' nanopore.fastq.gz
18 长度分布直方图
seqkit watch -L -f ReadLen hairpin.fa
19 平均质量直方图
seqkit watch -p 500 -O qhist.pdf -f MeanQual nanopore.fastq.gz
20 选取固定范围
seqkit range -r 200:300 nanopore.fastq.gz
21 移除重名序列
seqkit rmdup -n nanopore.fastq.gz -o clean.fa.gz
22 将小于 Q20 的替换为小写字母
seqtk seq -q 20 illumina_1.fastq.gz
写在最后:本片推文涉及的fastq.gz文件,在之前的推文都有介绍大家去下载。大家可以找任意一个来代入命令学习下fastq文件到底长什么样子。下面放之前推文的链接。