【直播】我的基因组（十一）:测序数据的比对
上一次直播中，我们对拿到手的测序数据进行了质控，测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了，这里我们先用参考基因组hg19，大家可以参照【直播】我的基因组（五）:测试数据及参考基因组的准备来下载参考基因组hg19，我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对，可以参考【直播】我的基因组（四）：计算资源的准备来安装，以及对hg19建立索引。
我们首先简单讲一下为什么要进行比对以及比对过程是怎样的？
可以看到我们到手的测序数据格式是fastq，每条reads都是150个碱基组成，如果只看这fastq，我们没办法得知它的意义，参考基因组那么大(人类约30亿个碱基)，这个reads在我们基因组的哪里呢？
简单解释一下，假设人类基因组是123456789，如果我们的测序得到的reads是123，那么我们很明显知道这条reads在基因组的第一个位置，跨越了3个长度，如果我们的reads是567，那么我们也可以看到它在基因组的第5个位置。如果我们看到了一个reads是235567，同样我们也很容易看到它在基因组第2位置，但是在第3个位置参考是4，它却是5，这里可能是测序错误，也可能是这个reads真的变异了！
但是我们的参考基因组远远没有那么简单，30亿个碱基的庞大数目，测序的一条reads也有150个碱基，仅仅用肉眼观察基本不可能判断出它到底在哪里。但并非一定观察不到，如果你有多的不可计的时间及精力的话，手工比对穷极一生来搞定一条reads的比对就很不容易了（当然肯定不会有人这么傻，这里只是说数据量真的很大而已）。然而在我们手上可是有8.9亿条reads，所以我们需要借助计算机来进行比对，现在比较流行的基因组比对工具是bwa和bowtie，它俩的算法不一样，但是我们不需要了解那么具体，只需要知道它可以把我们的fastq测序文件通过与参考基因组的比对生成sam格式（自行搜索了解该格式）的比对结果文件（如下），从sam文件中，我们可以看到每条reads在参考基因组的位置，这条reads是在哪一条染色体，又是在这条染色体的哪个位置就可以一目了然。
对于比对的结果，我们可以用IGV可视化查看，还可以手动查看每个基因的比对情况：
下面我简单讲一下代码
1，下载hg19基因组
cd ~/reference
mkdir -p genome/hg19 && cd genome/hg19
nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &
tar zvfx chromFa.tar.gz
cat *.fa > hg19.fa
首先要理解linux基础命令，在我们的服务器上面新建好目录，找到hg19的下载链接，用linux自带的wget下载，因为文件太大，所以我们用nohup放在后台下载。下载后是压缩文件 chromFa.tar.gz，在linux里面需要用tar zvfx 来解压tar.gz文件即可。解压开后是一个个文件，需要用cat合并！最终效果如下：
2，安装bwa软件
## Download and install BWA
cd ~/biosoft
mkdir bwa && cd bwa
#http://sourceforge.net/projects/bio-bwa/files/
wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2
tar xvfj bwa-0.7.15.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files
cd bwa-0.7.15
make
~/biosoft/bwa/bwa-0.7.15/bwa
我所有的软件都安装在自己的home目录下面的biosoft文件夹。同样，也是找的bwa的下载地址，然后解压，然后直接make即可。很多人的服务器会报错zlib.h缺少的问题，看我以前的教程：http://www.bio-info-trainee.com/518.html ,缺少什么你就安装什么，但是缺少的东西需要安装到系统环境变量，但是我的bwa是直接安装到自己的目录，所以我用全路径在调用该软件。如果你的这个命令~/biosoft/bwa/bwa-0.7.15/bwa 能够显示下面的help文档，说明你已经安装成功啦~
3，对hg19参考基因组用bwa构建索引
cd ~/reference
mkdir -p index/bwa && cd index/bwa
nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index -a bwtsw -p ~/reference/index/bwa/hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
代码很简单，就是新建好一个文件夹来存放我们的参考基因组的索引，我这里选择的是我的home目录下面的reference/index/bwa/ 文件夹，可以看到如下内容：
我还是用了nohup把这个命令挂在后台，防止掉线，因为要运行2个小时左右，我加上time命令可以看到运行时间，我用了bwa的index模式来索引参考基因在，具体bwa用法可以自己看文档，但是我们只需要学会索引及比对就好了。有点类似于window下面的软件有一个个菜单栏一样，需要自己的鼠标点击来实现一个个功能，在linux下面就是把命令准备好，然后运行。
4.把fastq文件比对到参考基因组
这个命令就一句话，但是里面的信息量非常大， 需要熟练掌握linux命令以及shell脚本的语法，但是解析起来也很简单，就是因为我们的fastq文件命名是有规律的，根据规律我构造出一个循环命令，里面的i这个变量会自动扩展成1，2，3，4，5，6依次来用bwa mem 模式来比对，因为是PE150测序，所以选择这个模式，-M就是选择我们上一步构建好的参考基因组，最后面的 1> 和2>是把软件运行结果输出来，分别是标准输出和标准错误输出，大家可以自行搜索。如果fastq文件的命名发生变化，这个shell脚本是运行不了的，需要临时构建，自己得掌握脚本编写，不然就一个个的比对，手动。
大家可以去看【直播】我的基因组（七）:从整体理解全基因组测序数据的变异位点，来了解这个命令的运行结果。
文：Jimmy、吃瓜群众
图文编辑：吃瓜群众