一、FASTA和FASTQ
1. FASTA
FASTA文件主要由两个部分构成：序列头信息（有时包括一些其它的描述信息）和具体的序列数据。头信息独占一行，以大于号（>）开头作为识别标记，其中除了记录该条序列的名字之外，有时候还会接上其它的信息。紧接的下一行是具体的序列内容，直到另一行碰到另一个大于号（>）开头的新序列或者文件末尾。下面给出一个FASTA文件的例子
>ENSMUSG00000020122|ENSMUST00000138518
CCCTCCTATCATGCTGTCAGTGTATCTCTAAATAGCACTCTCAACCCCCGTGAACTTGGT
TATTAAAAACATGCCCAAAGTCTGGGAGCCAGGGCTGCAGGGAAATACCACAGCCTCAGT
TCATCAAAACAGTTCATTGCCCAAAATGTTCTCAGCTGCAGCTTTCATGAGGTAACTCCA
GGGCCCACCTGTTCTCTGGT
>ENSMUSG00000020122|ENSMUST00000125984
GAGTCAGGTTGAAGCTGCCCTGAACACTACAGAGAAGAGAGGCCTTGGTGTCCTGTTGTC
TCCAGAACCCCAATATGTCTTGTGAAGGGCACACAACCCCTCAAAGGGGTGTCACTTCTT
CTGATCACTTTTGTTACTGTTTACTAACTGATCCTATGAATCACTGTGTCTTCTCAGAGG
CCGTGAACCACGTCTGCAAT
注意
- 除了序列内容之外，FASTA的头信息并没有被严格地限制。这个特点有时会带来很多麻烦的事情，比如有时我们会看到相同的序列被不同的人处理之后、甚至是在不同的网站上或者数据库中它们的头信息都不尽相同，比如以下的几种情况都是可能存在的。
>ENSMUSG00000020122|ENSMUST00000125984
> ENSMUSG00000020122|ENSMUST00000125984
>ENSMUSG00000020122|ENSMUST00000125984|epidermal growth factor receptor
>ENSMUSG00000020122|ENSMUST00000125984|Egfr
>ENSMUSG00000020122|ENSMUST00000125984|11|ENSFM00410000138465
这对于程序处理来说，凌乱的格式显然是不合适的。因此后来在业内也慢慢地有一些不成文的规则被大家所使用，那就是，用一个空格把头信息分为两个部分：第一部分是序列名字，它和大于号（>）紧接在一起；第二部分是注释信息，这个可以没有，就看具体需要。很多生信软件（如：BWA，samtools，bcftools，bedtools等）都是将第一个空格前面的内容认定为序列名字来进行操作的。
- FASTA由于是文本文件，它里面的内容是否有重复是无法自检的，在使用之前需要我们进行额外的检查。这个检查倒不用很复杂，只需检查序列名字是否有重复即可。但对于那些已经成为标准使用的参考序列来说，都有专门的团队进行维护，因此不会出现这种内容重复的情况，可以直接使用，但对于其它的一些序列来说，谨慎起见，最好进行检查。
2. FASTQ
上面所讲的FASTA文件，它所存的都是已经排列好的序列（如参考序列），FASTQ存的则是产生自测序仪的原始测序数据，它由测序的图像数据转换过来，也是文本文件，文件大小依照不同的测序量（或测序深度）而有很大差异，小的可能只有几M，大的则常常有几十G上百G，文件后缀通常都是.fastq，.fq或者.fq.gz（gz压缩），以下是它的一个例子：
@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码表示。
详细谈谈FASTQ质量值的计算方法
在测序仪进行测序的时候，会自动根据荧光信号的强弱给出一个参考的测序错误概率（error probility，P）根据定义来说，P值肯定是越小越好。我们怎么储存他们呢？直接储存成小数点？比如1%储存成0.01？这肯定是不高效的，因为1个碱基的信息，占用了至少4个字符。
所以科学家们的做法想了一个办法：
1.将P取log10之后再乘以-10，得到的结果为Q。
比如，P=1%，那么对应的Q=-10*log10（0.01）=20
2.把这个Q加上33或者64转成一个新的数值，称为Phred，最后把Phred对应的ASCII字符对应到这个碱基。
如Q=20，Phred = 20 + 33 = 53，对应的符号是”5”
这样就可以用1个符号与1个碱基一一对应，是不是很聪明？但一开始对于要加哪一个整数，并没有什么指导标准，这就导致了在刚开始的时候，不同的测序平台加的整数也不同，总的来说有以下3种质量体系，演变到现在也基本只剩下第一种（Phred33）了，如下表：
从表中可以看到下限有33和64两个值，我们把加33的的质量值体系称之为Phred33，加64的称之为Phred64（Solexa的除外，它叫Selexa64）。不过，现在一般都是使用Phred33这个体系。
参考:
https://zhuanlan.zhihu.com/p/20714540
二、FastQC、Multiqc查看质量
1. FastQC
FastQC是一款基于Java的软件，一般都是在linux环境下使用命令行运行，它可以快速多线程地对测序数据进行质量评估（Quality Control），其官网地址为：Babraham Bioinformatics
FastQC的下载和安装，和一般的Java软件没有什么区别，我们在这里就不做介绍了，在成功安装好以后，我们就在命令行模式下，输入fastqc就可以调用这个程序。
# 基本格式
# fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
# 主要是包括前面的各种选项和最后面的可以加入N个文件
# -o --outdir FastQC生成的报告文件的储存路径，生成的报告的文件名是根据输入来定的
# --extract 生成的报告默认会打包成1个压缩文件，使用这个参数是让程序不打包
# -t --threads 选择程序运行的线程数，每个线程会占用250MB内存，越多越快咯
# -c --contaminants 污染物选项，输入的是一个文件，格式是Name [Tab] Sequence，里面是可能的污染序列，如果有这个选项，FastQC会在计算时候评估污染的情况，并在统计的时候进行分析，一般用不到
# -a --adapters 也是输入一个文件，文件的格式Name [Tab] Sequence，储存的是测序的adpater序列信息，如果不输入，目前版本的FastQC就按照通用引物来评估序列时候有adapter的残留
# -q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况。
以我平时用的一个真实的例子：
fastqc -o ./tmp.result/fastQC/ -t 6 ./tmp.data/fastq/H1EScell-dnase-2014-GSE56869_20151208_SRR124
运行一段时间以后，就会出现报告：
H1EScell-dnase-2014-GSE56869_20151208_SRR1248176_1.fq_fastqc.html
H1EScell-dnase-2014-GSE56869_20151208_SRR1248176_1.fq_fastqc.zip
使用浏览器打开后缀是html的文件，就是图表化的fastqc报告：
2. Multiqc 一个高颜值的质量查看工具
fastqc是一款基于java的软件，能够对测序数据的质量进行评估。一个样本生成一个报告，当样本量过多时，逐一查看样本质量就稍显不方便，multiqc是一个基于Python的模块, 用于整合其它软件的报告的软件，能将fastqc生成的多个报告整合成一个报告的软件，这样能方便的查看所有测序数据的质量。
2.1 multiqc的安装：
conda install -c biocondamultiqc
2.1 现在用最简单的命令整合fastqc的报告：
(multiqc+fastqc结果报告存放路径+multiqc报告输出路径）
multiqc /data/home/chj/fastqc_result -o/data/home/chj/multiqc_result
命令执行完毕会生成1个html报告，直接网页打开就可以查看和一个multiqc_data的文件夹，其中包含一些数据基本的统计信息和日志文档。
参考：
https://zhuanlan.zhihu.com/p/20731723
https://www.jianshu.com/p/85da4dcc6020
三、质量控制与过滤
我们已经知道现在的NGS测序，以illumina为首基本都是运用边合成边测序的技术。碱基的合成依靠的是化学反应，这使得碱基链可以不断地从5'端一直往3'端合成并延伸下去。但在这个合成的过程中随着合成链的增长，DNA聚合酶的效率会不断下降，特异性也开始变差，这就会带来一个问题——越到后面碱基合成的错误率就会越高，这也是为何当前NGS测序读长普遍偏短的一个原因。
1. Trimmomatic切除测序接头序列和read的低质量序列
目前也已有很多工具用来切除接头序列和低质量碱基，比如SOAPnuke、cutadapt、untrimmed等不下十个，但这其中比较方便好用的是Trimmomatic（也是一个java程序）、sickle和seqtk。
Trimmomatic的好处在于，它不但可以用来切除illumina测序平台的接头序列，还可以去除由我们自己指定的特定接头序列，而且同时也能够过滤read末尾的低质量序列，sickle和seqtk只能去除低质量碱基。具体的原理就是通过滑动一定长度的窗口，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除，注！意！不是挖掉read中的这部分低质量序列，而是像切菜一样，直接从低质量区域开始把这条read后面的所有其它碱基全！部！剁！掉！否则就是在人为改变实际的基因组序列情况。
1.1 首先是安装Trimmomatic。我们可以到它的官网上获取最新的版本，下载打包好的binary即可]
$ java -jar trimmomatic-0.36.jar
同个目录下还有一个名为adapters的文件夹，这个文件夹中的内容对于我们去除接头序列来说非常重要。其中默认存放的是illumina测序平台的接头序列（fasta格式），在实际的使用过程中，如果需要去除接头，我们需要明确指定对应的序列作为输入参数。
1.2 Trimmomatic有两种运行模式：PE和SE。顾名思义，PE就是对应Pair End测序的，SE则是对应Single End测序的。
PE模式，HiSeq PE测序：
$ java -jar /path/Trimmomatic/trimmomatic-0.36.jar PE -phred33 -trimlog logfile reads_1.fq.gz reads_2.fq.gz out.read_1.fq.gz out.trim.read_1.fq.gz out.read_2.fq.gz out.trim.read_2.fq.gz ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
SE模式，HiSeq SE测序：
java -jar /path/Trimmomatic/trimmomatic-0.36.jar SE -phred33 -trimlog se.logfile raw_data/untreated.fq out.untreated.fq.gz ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:50
1.3 参数说明
ILLUMINACLIP：接头序列切除参数。LLUMINACLIP:TruSeq3-PE.fa:2:30:10（省掉了路径）意思分别是：TruSeq3-PE.fa是接头序列，2是比对时接头序列时所允许的最大错配数；30指的是要求PE的两条read同时和PE的adapter序列比对，匹配度加起来超30%，那么就认为这对PE的read含有adapter，并在对应的位置需要进行切除。10和前面的30不同，它指的是，反正只要这条read的某部分和adpater序列有超过10%的匹配率，那么就代表含有adapter了，需要进行去除；
【注】测序的时候一般只会测到一部分的adapter，因此read和adaper对比的时候肯定是不需要要求百分百匹配率的，上述30%和10%其实是比较推荐的值。
SLIDINGWINDOW：滑动窗口长度的参数，SLIDINGWINDOW:5:20代表窗口长度为5，窗口中的平均质量值至少为20，否则会开始切除；
LEADING：规定read开头的碱基是否要被切除的质量阈值；
TRAILING：规定read末尾的碱基是否要被切除的质量阈值；
MINLEN：规定read被切除后至少需要保留的长度，如果低于该长度，会被丢掉。
2. 使用fastp进行数据质控
2.1简单介绍
fastp是一款较新的数据质控软件，并且运行的速度要快于基于java的Trimmomatic，fastp软件会生成HTML格式的报告，而且该报告中没有任何一张静态图片，所有的图表都是使用JavaScript动态绘制，非常具有交互性。想要看一下样板报告的，可以去以下链接：
http://opengene.org/fastp/fastp.html
而且软件的开发者还充分考虑到了各种自动化分析的需求，不但生成了人可读的HTML报告，还生成了程序可读性非常强的JSON结果，该JSON报告中的数据包含了HTML报告100%的信息，而且该JSON文件的格式还是特殊定制的，不但程序读得爽，你用任何一款文本编辑器打开，一眼过去也会看得明明白白。想要看一下JSON结果长什么样的，可以去以下链接：
http://opengene.org/fastp/fastp.json
2.2 具体用法
fastp -i 20-m-1_FRAS202147686-1r_1.fq.gz -I 20-m-1_FRAS202147686-1r_2.fq.gz -o 20m1_good_1.fq -O 20m1_good_2.fq -q 30 -l 100 -n 10 --detect_adapter_for_pe &
usage: fastp -i <in1> -o <out1> [-I <in1> -O <out2>] [options...]
options:
# I/O options 即输入输出文件设置
-i, --in1 read1 input file name (string)
-o, --out1 read1 output file name (string [=])
-I, --in2 read2 input file name (string [=])
-O, --out2 read2 output file name (string [=])
# length filtering options 根据序列长度来过滤序列
-L, --disable_length_filtering length filtering is enabled by default. If this option is specified, length filtering is disabled
-l, --length_required reads shorter than length_required will be discarded, default is 15. (int [=15])
# quality filtering options 根据碱基质量来过滤序列
-Q, --disable_quality_filtering quality filtering is enabled by default. If this option is specified, quality filtering is disabled
-q, --qualified_quality_phred the quality value that a base is qualified. Default 15 means phred quality >=Q15 is qualified. (int [=15])
-u, --unqualified_percent_limit how many percents of bases are allowed to be unqualified (0~100). Default 40 means 40% (int [=40])
-n, --n_base_limit if one read's number of N base is >n_base_limit, then this read/pair is discarded. Default is 5 (int [=5])
# help
-?, --help print this message
参考：
https://www.jianshu.com/p/6f492058da5b
网友评论