【直播】我的基因组(十一):测序数据的比对

上一次直播中,我们对拿到手的测序数据进行了质控,测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了,这里我们先用参考基因组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

rm chr*.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文件比对到参考基因组

for i in $(seq 1 6) ;do (nohup ~/biosoft/bwa/bwa-0.7.15/bwa  mem -t 5 -M ~/reference/index/bwa/hg19  KPGP-00001_L${i}_R1.fq.gz KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log &);done

这个命令就一句话,但是里面的信息量非常大, 需要熟练掌握linux命令以及shell脚本的语法,但是解析起来也很简单,就是因为我们的fastq文件命名是有规律的,根据规律我构造出一个循环命令,里面的i这个变量会自动扩展成1,2,3,4,5,6依次来用bwa  mem 模式来比对,因为是PE150测序,所以选择这个模式,-M就是选择我们上一步构建好的参考基因组,最后面的 1> 和2>是把软件运行结果输出来,分别是标准输出和标准错误输出,大家可以自行搜索。如果fastq文件的命名发生变化,这个shell脚本是运行不了的,需要临时构建,自己得掌握脚本编写,不然就一个个的比对,手动。

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

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

1

Comments are closed.