基于二代测序的序列比对的文件格式和含义

二代测序又叫做NGS（Next-generation-sequencing），是将生物体基因组打断后，用Illumina或Ion Torrent平台的测序仪进行碱基读取的技术。目前华大基因也已推出比较新型的二代测序仪，比如MGISEQ-2000， MGISEQ-T7，这些读取的DNA序列片段长度普遍在150-300bp之间。

不同生物体的基因组大小也不一样，含有的染色体数目也不同，这里以人的基因组为例。

可以从一个UCSC的网站下载   http://hgdownload.cse.ucsc.edu/downloads.html#human  （如下图，在Full data set中）

这里面有两个主要的人类基因组版本，hg19和hg38，区别主要在于基因组的精密度不一样，hg38采用了更多的人类DNA样本，拼接的更完整一些。不过现在的科研和基因检测两种版本都在用。

hg19的各个染色体的大小如下：

chr1 249250621

chr2 243199373

chr3 198022430

chr4 191154276

chr5 180915260

chr6 171115067

chr7 159138663

chr8 146364022

chr9 141213431

chr10 135534747

chr11 135006516

chr12 133851895

chr13 115169878

chr14 107349540

chr15 102531392

chr16 90354753

chr17 81195210

chr18 78077248

chr19 59128983

chr20 63025520

chr21 48129895

chr22 51304566

chrX 155270560

chrY 59373566

chrM 16569

采用的计算脚本如下：

#!/usr/bin/env python
# -*- coding : UTF-8 -*-

from Bio import SeqIO
import sys

fastafile = sys.argv[1] 

for s in SeqIO.parse(fastafile, "fasta"):
	if len(s.id) < 10:
		print("{0}\t{1}".format(s.id, len(s.seq)))

序列比对这一步主要就是将NGS得到的 150bp或300bp的reads与基因组的序列进行字符串比较，找出这些reads在基因组上的比对位置，以及测到的与基因组某些位置不一样的碱基信息（也就是基因检测领域常说的突变，其中包括SNP/INDEL，CNV和SV），分别叫做单核苷酸多态性，小片段的插入和缺失，拷贝数变异和结构变异。

下面讲一下NGS测序仪得到的序列文件格式，大多数测序仪生成的序列格式为FASTQ， 也有一些是BAM格式的。

这里主要讲一下FASTQ格式（如图）：

一个fastq文件中可能有几百万条序列，以上这张图片中只展示了三条序列，每条序列占了四行。

第一条序列的ID是FP1LH:08681:03840

序列的碱基排列就是第二行的内容

第三行是占位符，没有实际含义

第四行是这条序列的各个碱基在测序的时候测到的质量值的ASCII码转换形式

第一条序列在hg19上的比对位置如下 （这里我才用UCSC 的blat工具比对得出）：

blat的网址为  http://genome.ucsc.edu/cgi-bin/hgBlat

可以看出这条reads 比对到了1号染色体， 第903017到903144个base pair这一段上，而且比对到了DNA 的负链上。

如果是直接用fastq文件去与hg19的fasta进行比对，就不能用blat这个工具了，一般要用bwa或者bowtie这样的专业二代测序比对软件(在Linux操作系统下可用，需要把软件的路径添加到环境变量)。

bwa的使用方法如下(比对前需用bwa index命令对hg19.fasta构建index，耗时约2h)：

bwa mem -p -t 8 hg19.fasta test.fastq > output.sam 2>./log/align_bwa.log

bowtie的使用方法如下

(需要用bowtie-build对hg19.fasta构建比对的index)

bowtie -p 8 -l 5 --chunkmbs 512 -v 3 -k 50 hg19.fasta test.fastq -S output.sam 1>./log/align_bowtie.log 2>&1

比对的结果都输出到SAM格式中，SAM格式参考 http://samtools.sourceforge.net/SAM1.pdf

如果是单端测序数据，每条reads的比对信息只占一行，如果是双端测序数据，每对reads比对信息占两行。

在上图中可以看到，FP1LH:08681:03840 这条reads 比对的flag值为16，flag的值的含义可以参考这里（https://broadinstitute.github.io/picard/explain-flags.html）,  在SAM flag 输入框里输入一个flag值，点击Explain，可以得到该值的含义。

比对位置在chr1的903017， 比对质量值为60，这个值是很高的。reads长度为128bp，而且全都Match。

其它各列的含义可以从互联网中查到，这里就不做解释了。