上一次直播中，我们对拿到手的测序数据进行了质控，测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了，这里我们先用参考基因组hg19，大家可以参照【直播】我的基因组（五）:测试数据及参考基因组的准备来下载参考基因组hg19，我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对，可以参考【直播】我的基因组（四）：计算资源的准备来安装，以及对hg19建立索引。

可以看到我们到手的测序数据格式是fastq，每条reads都是150个碱基组成，如果只看这fastq，我们没办法得知它的意义，参考基因组那么大(人类约30亿个碱基)，这个reads在我们基因组的哪里呢？

但是我们的参考基因组远远没有那么简单，30亿个碱基的庞大数目，测序的一条reads也有150个碱基，仅仅用肉眼观察基本不可能判断出它到底在哪里。但并非一定观察不到，如果你有多的不可计的时间及精力的话，手工比对穷极一生来搞定一条reads的比对就很不容易了（当然肯定不会有人这么傻，这里只是说数据量真的很大而已）。然而在我们手上可是有8.9亿条reads，所以我们需要借助计算机来进行比对，现在比较流行的基因组比对工具是bwa和bowtie，它俩的算法不一样，但是我们不需要了解那么具体，只需要知道它可以把我们的fastq测序文件通过与参考基因组的比对生成sam格式（自行搜索了解该格式）的比对结果文件（如下），从sam文件中，我们可以看到每条reads在参考基因组的位置，这条reads是在哪一条染色体，又是在这条染色体的哪个位置就可以一目了然。

对于比对的结果，我们可以用IGV可视化查看，还可以手动查看每个基因的比对情况：

下面我简单讲一下代码

1，下载hg19基因组

首先要理解linux基础命令，在我们的服务器上面新建好目录，找到hg19的下载链接，用linux自带的wget下载，因为文件太大，所以我们用nohup放在后台下载。下载后是压缩文件 chromFa.tar.gz，在linux里面需要用tar zvfx 来解压tar.gz文件即可。解压开后是一个个文件，需要用cat合并！最终效果如下：

2，安装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构建索引

代码很简单，就是新建好一个文件夹来存放我们的参考基因组的索引，我这里选择的是我的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脚本是运行不了的，需要临时构建，自己得掌握脚本编写，不然就一个个的比对，手动。

大家可以去看【直播】我的基因组（七）:从整体理解全基因组测序数据的变异位点，来了解这个命令的运行结果。

请扫描以下二维码关注我们，获取直播系列的所有帖子！