068-短序列比对了解一下?

刘小泽写于18.12.17

组学分析离不开比对环节,可能流程已经会了,但这不是根本。我们需要再深入学习一下关于短序列比对的一些知识,这样以后再分析的时候,就更能明白我们自己在做些什么

短序列是什么?

一代测序一般都测的比较长(1000bp左右)并且准确率很高,能到99.9%,但最大的缺点就是:这么好的技术却不能大批量应用,只能一条一条测,耗时耗钱。从2005年左右,二代测序的短读长、高通量开始迅速发展,产生的reads和一代有着显著的不同:序列长度明显降下来了(50-300bp),一次测序产量明显上去了(一次上亿条reads)。这里不讨论测序仪怎么做到的输出短序列,只需要知道二代测序带给我们的结果一是短,二是多就行了。

有了输出的短序列,然后呢?

它们的用处一般有两个:一个是比对到参考基因组(也就是重测序 re-sequencing流程),另一个就是相互比对进行叠加(也就是组装的过程,reads=》contig=》scaffold等)。这些工作都是靠比对工具完成的,我们称之为:“short read aligners” OR “short read mappers”

比对工具怎么完成的?

首先,让我们感谢那些**短序列比对工具(以下简称“工具”)**的作者,生信界有了他们才有了不同组合的流程,才推动了不同的组学研究。他们设计的一个好的工具可以在1秒内比对超过一万条序列到人类基因组的30亿bp碱基上

比对一般需要寻找exact matchesmissing locations 以及partial and alternative matches/overlaps ,但是这几个任务每个软件都有自己的擅长与不足。另外,根据数据类型(domain-specific)表现也不同

当然,工具也有自己的劣势,大部分的工具只能根据阈值去寻找和目标区域相似的比对结果,高相似性优先,因此有部分会因阈值设定被抛弃;当序列过长或过短时,就会受到限制,例如研究small RNA这种低于30nt的序列,可选的比对工具就是少之又少,因为对比对算法的精确度要求更高

alignment vs mapping

说起比对,我们一般都会想起来这两个词,但是严格来说,它们还是有区别的:mapping强调的是不要求比对的准确度,能有比对就可以,是region层面上的,它的目的就是看看基因组上有没有一块区域可以安放这条序列;而alignment强调精确比对,是碱基层面上的

李恒关于这二者做过自己的解释:

Mapping

  • A mapping is a region where a read sequence is placed.
  • A mapping is regarded to be correct if it overlaps the true region.

Alignment

  • An alignment is the detailed placement of each base in a read.
  • An alignment is regarded to be correct if each base is placed correctly.

例如:如果研究SNPs和基因组内变异,应该就是用alignment的工具;如果研究RNA-seq就是用mapping的工具,因为RNA-seq一般我们就想知道这个序列从哪里来的就好了,没必要知道比对上多少

这么多工具,如何选择?

其实工具很多,能被大众接受的很少,因此软件开发之路十分艰难,其实大部分的开发者并没有得到关注。

这里用的最多的当属bwabowtie2 ,选择的话一般就是根据喜好和信仰充值了。我们知道bwa是曾经的博德研究所的李恒开发的,国人的肯定要支持;另外华大也开发了一个novoalign 喜欢BGI的也可以关注。但是不管怎样,一个比对项目,可以多用几种工具分析分析,增加自己结果的可信度

另外,在选择工具的时候,我们要根据自己的实际需求,自己做一个心理暗示(mental checklist):

  • 比对算法:global, local or semi-local
  • 工具遇到INDEL怎么处理?
  • 工具可以跨大片段区域比对吗?
  • 工具可以进行根据阈值过滤吗?阈值设置?
  • 工具可以找到嵌合体比对吗?

选择一个最常用的bwa研究下

bwa的背景

bwa是Burrows-Wheelers Aligner的缩写,应该是目前使用最为广泛的,后期会搭配同作者的samtools使用

当时李恒在发表这个bwa算法的时候,那时已经好多人在用了,还被拒稿,为此他写过一封信:

img

生信有一个特点,就是:更新速度快,可能有的软件还不知道名字就已经被淘汰了。身处前端的bwa也是不断的更替版本,之前的算法是 bwa aln,现在是bwa mem【这里的mem指的是Maximally Exact Matches】。因此如果采用不同的算法,虽然都是bwa,但结果还是有差别。这就是有时让人纠结的地方,因为不清楚作者软件的版本更替,于是,最好的办法就是用熟练github,自己去找(当然,我也在努力的道路上)

bwa准备工作

一般来讲,使用bwa比对之前需要根据参考基因组构建索引index ,目的就是快速搜索基因组。一般参考基因组放一个目录,然后构建的index放另一个目录,随时取用。

bwa小练习

将埃博拉病毒测序的reads比对到它的参考基因组上(命名为1976.fa)

先构建索引
mkdir -p ~/tmp/bwa_test/ref

# 下载基因组(19Kb)[安装entrez-direct]
conda install -c bioconda entrez-direct
efetch -db=nuccore -format=fasta -id=AF086833 > ~/tmp/bwa_test/ref/1976.fa

# 构建索引
ref=~/tmp/bwa_test/ref/1976.fa
bwa index $ref

# ls ~/tmp/bwa_test/ref
然后下载测试数据
mkdir -p ~/tmp/bwa_test/raw && cd ~/tmp/bwa_test/raw

# 获取全部的埃博拉病毒项目的测序数据
esearch -db sra -query PRJNA257197 | efetch -format runinfo > runinfo.csv

# 挑选SRR1972739下载
cat runinfo.csv| grep SRR1972739 | cut -d "," -f 10 | xargs -n1 wget {}

# 解压成fq文件 [为了测试,选取前1万条reads]
fastq-dump -X10000 --split-files SRR1972739
比对

比对默认得到sam格式(Sequence Alignment Map),其中包括了所有的比对信息,之后再也不用看fastq了

mkdir -p ~/tmp/bwa_test/align && cd ~/tmp/bwa_test/align

r1=~/tmp/bwa_test/raw/SRR1972739_1.fastq
r2=~/tmp/bwa_test/raw/SRR1972739_2.fastq
ref=~/tmp/bwa_test/ref/1976.fa

baw mem $ref $r1 $r2 > bwa.sam

至于Bowtie2,也有自己的用法,不过是将索引换成了bowtie2-build , 然后比对的时候bowtie2 -x $ref -1 $r1 -2 $r2 > bowtie.sam

运行结束后,会给个报告

软件的比较?

如果只考虑时间的话,bwa要快一些,当然这也是有原因的:bowtie2产生的结果中包含有更多的信息,例如,同时看下它们的tag部分:

cat bwa.sam | cut -f 12-20 | head

cat bowtie.sam | cut -f 12-20 | head

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related