2020-01-11 了解FASTQ格式并处理FASTQ文件

FASTQ文件格式是测序仪展示数据的标准格式，可以看成FASTA文件的变种（FASTA+Q），因为其包含了对序列中每个碱基的Qualify Measurement。（如：碱基A出错的可能性是1/1000）

FASTQ格式详述

FASTQ格式包括4个部分，每个部分1行，格式同FASTA相似，但缺陷也更多：

1. 类似FASTA的头部，以@而非>起始，后跟ID和描述文本
2. 测定的序列，通常为1行，但有时也会换行，最后以+指示下一部分
3. 由+表示（后面有时会跟着和第一部分相同的id和header）
4. 编码第2部分测定序列的质量值，长度必须同第2部分相同，换行方式也要同第2部分相同
   
   第4部分看着有点奇怪，其实是通过转码将两位数字的Phred Score转换为1个字符的Quality Score

第一行为FASTQ quality codes
第二行为Quality Score (Q)/Phred Score (P)

Sanger(+33)格式

错误率公式：Error=10ˆ(-P/10)

编码为I，P=40，错误率为10^(-40/10)=0.01%

以前还用过一种老的+64格式的FASTQ编码：

如果你的FASTQ文件中Quality Score为小写，代表你的FASTQ是老版本
可以使用seqtk命令进行转换
seqtk seq -Q64 input.fq > output.fa

FASTQ文件header中的额外信息

参考维基百科fastq format词条
EAS139：仪器，
136：run编号，
FC706VJ：流通池编号，
2：泳道，
2104：tile编号，
(15343,197393)：xy坐标，
1: the member of a pair, 1 or 2 (paired-end or mate-pair reads only)
Y: Y if the read is filtered, N otherwise
18: 0 when none of the control bits are on, otherwise it is an even number
ATCACG指索引序列
此信息针对特定仪器/供应商，可能会随仪器的不同版本或版本而变化。但对于回答可能影响仪器的质量控制或系统性技术问题是很有用的。

通过命令行转换FASTQ质量编码

略，在书p185

进阶FASTQ处理

介绍使用SeqKit操作FASTA/Q，SeqKit支持FASTA/Q文件和FASTA/Q的压缩格式：

seqkit命令

seqkit的学习笔记可参考

《fasta/fq文件处理万能工具——Seqkit学习记录》https://www.jianshu.com/p/f0e65738b7c7

还需要使用csvtk（CSV/TSV）：
下载
https://github.com/shenwei356/csvtk/releases
安装csvtk的方式：
https://bioinf.shenwei.me/csvtk/#installation

#作者目前尚不推荐使用conda进行安装
curl http://app.shenwei.me/data/csvtk/csvtk_linux_amd64.tar.gz > csvtk_linux_amd64.tar.gz
tar -zxvf csvtk_linux_amd64.tar.gz
sudo cp csvtk bin

1. 示例数据

使用一个FASTQ文件（Sample Reads，1M）和两个FASTA文件（NCBI RefSeq数据库中的病毒DNA和蛋白质序列，60M+40M）。

#wget -nc指若要覆盖已有文件则不要下载
wget -nc http://data.biostarhandbook.com/reads/duplicated-reads.fq.gz
wget -nc ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.2.1.genomic.fna.gz
wget -nc ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.2.protein.faa.gz
# fna文件和faa文件
# *.fna = FASTA Nucleic Acid file
# *.faa = FASTA Amino Acid file

2. FASTQ文件概览

seqkit stat *.gz

3. 得到GC含量

seqkit fx2tab将FASTA/Q转换为3列表格形式（1:名称/ID，2:序列，3:质量），还可以在新列中提供各种信息，包括序列长度、GC内容/GC skew、alphabet
seqkit fx2tab -n -i -g viral.2.1.genomic.fna.gz | head
-n, --name only print names (no sequences and qualities)
-i, --only-id print ID instead of full head
-g, --gc print GC content
可以尝试分别去掉这些参数

4. 得到感兴趣碱基的含量

假设想要获得A、C、A+C的含量
seqkit fx2tab -H -n -i -B a -B c -B ac viral.2.1.genomic.fna.gz | head -5
-H, --header-line输出header line
-B, --base-content strings 输出碱基含量，忽略大小写，支持N等alphabet

5. 抽出部分序列和name/list文件

seqkit sample 根据数量或比例对序列进行抽样
-n, --number int 根据数量抽样（匹配可能不精确）
-p, --proportion float 根据比例抽样

seqkit seq seq命令对序列进行操作， (包括revserse, complement, extract ID等...)
-n, --name 只输出名字
-i, --only-id 只输出id而不输出full head

seqkit sample --proportion 0.001 duplicated-reads.fq.gz | seqkit seq --name --only-id > id.txt
head id.txt #id.txt中的记录要用于seqkit grep的检索，格式为每行1个记录

seqkit grep 根据ID/名称/序列/序列motif检索序列，允许错配
-f, --pattern-file string pattern file (one record per line)

seqkit grep --pattern-file id.txt duplicated-reads.fq.gz > duplicated-reads.subset.fq.gz
#用id.txt中的pattern在duplicated-reads.fq.gz中进行匹配

6. 找到包含简并碱基（degenerate base）的序列并定位

seqkit fx2tab --name --only-id --alphabet  viral.1.1.genomic.fna.gz | csvtk --no-header-row --tabs grep --fields 4 --use-regexp --ignore-case --pattern "[^ACGT]"

代码的长参数版本

seqkit fx2tab -n -i -a viral.2.1.genomic.fna.gz | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]"

seqkit fx2tab将FASTA / Q转换为表格格式，并可以在新列中输出序列字母:
-a, --alphabet 输出alphabet
然后可以使用文本搜索工具来过滤表。
csvtk下的参数
-H, --no-header-row 输入的CSV文件没有header行
-t, --tabsinput CSV file由tabs分隔。覆盖-d和-D

csvtk grep下的参数
-f, --fields string comma separated key fields, column name or index. e.g. -f 1-3 or -f id,id2 or -F -f "group" (default "1") 这个参数没看懂
-r, --use-regexp 给出的pattern为正则表达式
-i, --ignore-case 忽略大小写
-p, --pattern strings query pattern (multiple values supported) 这个参数没看懂
pattern为"[^ACGT]" 是一个正则表达式 这个正则表达式没看懂

参考博文：《csvtk——CSV_TSV文本处理万能工具》
https://www.jianshu.com/p/3aa6231ed5c3

保存这些包含简并碱基的序列的ID

seqkit fx2tab -n -i -a viral.2.1.genomic.fna.gz | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]" | csvtk -H -t cut -f 1 > id2.txt

使用seqkit grep命令去除

seqkit grep --pattern-file id2.txt --invert-match viral.1.1.genomic.fna.gz > clean.fa
#或
seqkit grep -P id2.txt --invert-match viral.1.1.genomic.fna.gz > clean.fa

csvtk cut下的参数
-f, --fields string 仅保留fields中的列，例： -f 1,2即保留第1、2列

7. 去除包含重复序列的FASTA/Q记录

8. 定位FASTA/Q序列中的模体/subsequence/酶消化位点

9. 根据序列长度排列FASTA/Q序列

10. 根据header信息spilt FASTA序列

11. 使用文本文件中的character strings在FASTA header中搜索和替换？

12. 从两个paired end reads文件中提取paired reads？

13. 连接两个FASTA序列