近日来苏州阴雨连绵，导致人不是很爽利。实验也并不顺利，不懂的东西还太多，要学的东西还太多。不过还好身边有开心的事，有可爱的人——周四看明白马尔科夫链的矩阵形式很开心，前天有人劝我每天都要开心很开心，昨晚出去吃饭很开心，今天电梯终于好了也很开心。总而言之，还是开心居多吧。
说了这么多闲话，下面进入正题。在我们终于拿到了fastq文件之后，又该做些什么呢？
1、fastq文件简介
1.1、格式简介
fastq格式是一种包含质量值的序列文件，其中的q为quality,一般用来存储原始测序数据，扩展名一般为fastq或者fq。目前illumina测序，BGISEQ，Ion Torrent，pacbio，nanopore都以fastq格式存储测序数据，其中illumina，BGISEQ一般是双末端测序，一般是一对文件，命名为_R1.fq.gz与_R2.fq.gz。下面是fastq格式常见的序列格式。
代码语言：javascript复制@FCD056DACXX:3:1101:2163:1959#TCGCCGTG/1
TCCGATAACGCTCAACCAGAGGGCTGCCAGCTCCGATCGGCAGTTGCAACCCATTGGCCGTCTGAGCCAGCAACCCCGGA
gggiiiiiiiiiiiiiiiiiiiiiiiiiigggggeeecccccc^bcbcccccccbccccc]aaccbbccc^R^^acccc_
@FCD056DACXX:3:1101:2194:1984#TCGCCGTG/1
AGACGACGACTTCGTTTCCCGCCGCGAGTTGCGCCATGATCGCGGTGTGCAGATTCGTTACGCCCTGGGCCACGGAGACG
gggiihiiiiiiiihiiiiiiiiiigeccccccccccccccccccaccccdcccccccccccacc_accccccccccV^^第一行：以‘@’开头，是这一条read的名字，这个字符串是根据测序时的状态信息转换过来的，中间不会有空格，它是每一条read的唯一标识符，同一份FASTQ文件中不会重复出现，甚至不同的FASTQ文件里也不会有重复；
第二行：测序read的序列，由A，C，G，T和N这五种字母构成，这也是我们真正关心的DNA序列，N代表的是测序时那些无法被识别出来的碱基；
第三行：以‘ ’开头，在旧版的FASTQ文件中会直接重复第一行的信息，但现在一般什么也不加（节省存储空间）；
第四行：测序read的质量值，这个和第二行的碱基信息一样重要，它描述的是每个测序碱基的可靠程度，用ASCII码表示。
1.2、质量值
上面提到fastq格式中的q代表质量值，因此fastq格式中质量值具有重要的作用，在很多的分析中会用到这个质量值，例如数据质控，数据过滤，序列拼接，短序列比对，变异检测中都要用到这个质量值。这个质量值是基于phred质量值体系，但是由于单个碱基无法与两位的质量值相匹配，例如A碱基对应的质量值为40，一个A字符对应两个字符40，因此需要将原始质量值加上33或者64，在转换为对应的ASCII码值，为何加33，因为33以下ASCII码无法用键盘字符表示出来。illumina测序1.8版本以上加33，以下加64。
质量值的计算方式是Q=-log10(错误率)*10，也就是错误率1%，Q为20。
1.3、完整性校验
完整性检验主要是为了保证文件在传输过程中保持完整，没有丢失内容，一般采用md5校验方式，目前测序公司给定的测序数据都带有md5文件，这样文件就是用来校验数据完整性的。可以使用md5sum -c命令检测这个文件，如果返回OK，说明文件完整。
代码语言：javascript复制md5sum -c SRR11178348_1.fastq.md5
md5sum -c SRR11178348_2.fastq.md5fastq文件格式简介部分来自以下文章，在此致谢。如有进一步需求，请移步原文。 王通，公众号：基因学苑fastq格式文件处理大全（一）
2、原始文件质控
2.1、单个文件质控
代码语言：javascript复制mkdir 3.fastq_qc
cd ./3.fastq_qc/
ln -s ~/path/to/fastq/*.fastq.gz ./
fastqc -t 16 -o ./ ./*.fastq.gz会得到一个zip文件和一个html报告
2.2、html报告解读
可以看到文件名、reads总数、读长是PE150，GC含量53%等信息。正常人类基因组的GC含量应该是43%左右，但是由于建库过程中不可避免的扩增偏倚等原因，最终的GC含量大概是50%左右。
测序质量尚可，没有小于20的，但是末端有一部分小于30，可以去除也可以不去除。
这里有点问题，虽然总体GC含量基本正常，但是在理论分布的右侧还有一个GC峰，可能是rRNA，因为这次的示例数据用的是去除rRNA建库的，而这种建库方式一般很难完全将rRNA除去
3、合并fastqc结果
3.1、multiqc合并qc结果
代码语言：javascript复制conda install -c bioconda multiqc
multiqc ./*zip -o ./3.2、查看html报告
可以看到基本所有的数据都有明显的GC双峰。此外，虽然对于RNA-seq来说，重复序列不可避免，但是高达80%的重复片段肯定也是有问题的，需要我们进一步对数据进行校正。
