搜索
查看: 4105|回复: 3

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

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-2 09:30:48 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-2-2 09:33 编辑

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

上一次直播中,我们对拿到手的测序数据进行了质控,测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了,这里我们先用参考基因组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基因组
[Perl] 纯文本查看 复制代码
cd ~/reference
mkdir -p genome/hg19  && cd genome/hg19
nohup wget [url=http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz]http://hgdownload.cse.ucsc.edu/g ... Zips/chromFa.tar.gz[/url] &
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软件

[Perl] 纯文本查看 复制代码
[mw_shl_code=perl,true]## Download and install BWA
cd ~/biosoft
mkdir bwa &&  cd bwa
#[url=http://sourceforge.net/projects/bio-bwa/files/]http://sourceforge.net/projects/bio-bwa/files/[/url]
wget [url=https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2]https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2[/url]
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构建索引
[Perl] 纯文本查看 复制代码
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文件比对到参考基因组

[Perl] 纯文本查看 复制代码
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脚本是运行不了的,需要临时构建,自己得掌握脚本编写,不然就一个个的比对,手动。


文:Jimmy、吃瓜群众

图文编辑:吃瓜群众





本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:小白生信学习记(二):服务器及其使用介绍
下一篇:【直播】我的基因组(十二):先粗略看看几个基因吧
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

15

帖子

88

积分

注册会员

Rank: 2

积分
88
发表于 2017-2-18 23:08:12 | 显示全部楼层
前面讲到会用BWA,BOWTIE2,HISAT做mapping,期待能讲到三者mapping的variant calling结果的区别及建议。
回复 支持 反对

使用道具 举报

0

主题

1

帖子

14

积分

新手上路

Rank: 1

积分
14
发表于 2017-2-25 21:56:30 | 显示全部楼层
Jimmy,不但给代码,还详细解释了代码,用心良苦,真心感谢~
回复 支持 反对

使用道具 举报

3

主题

11

帖子

258

积分

中级会员

Rank: 3Rank: 3

积分
258
发表于 2018-8-14 15:57:41 | 显示全部楼层
请问bwa mem 的 -M参数是什么用途?
-M        Mark shorter split hits as secondary (for Picard compatibility).
这个不太理解
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-20 07:22 , Processed in 0.031933 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.