二代测序call indel 总结

 

自己也做一代测序，二代测序数据分析3年有余，但是从来没有认真整理过；

看到这边博文写的不错，就转载过来了，先对做技术总结的的人致敬；

原文：http://blog.sina.com.cn/s/blog_6d25c0920101d3fn.html

Change logs:

Version 1.0.0 2012/12/08: The initial version.

 

Indel Calling相比于SNP Calling的难度要大一些，因为由于这种插入-缺失的存在，本身就很容易干扰排序，这种干扰会导致Indel周围出现很多假阳性的SNP，而且会影响Indel本身的准确性。理论上来说，检测Indel的最好方式就是做de novo assembly，然后比较de novo得到的基因组与原来的基因组，不过实际上de novo assembly的难度更大OTL

 

Paired-end测序为寻找较长片段的Indel提供了非常有用的信息，但是如何准确的利用这些信息也是目前这一块的一大难点。

 

下面就简单的介绍一些现有的Call Indel的软件及使用流程，注意一下这边介绍的几个基本上都是针对Illumina的paired-end数据的。

 

 

(1) Samtools mpileup

http://samtools.sourceforge.net/

 

Samtools里面的mpileup既可以Call SNP，也能Call Indel，

 

           samtools mpileup -DSugf ref.fasta 
                sample.bam | 
                bcftools view -Ncvg - 
                > sample.samtools.vcf

 

这边有个-N的参数表示跳过reference里面碱基是N的位点。

 

注意这边其实用的是一个管道，因为mpileup直接产生的结果不是最后vcf文件里看到的结果，而是一个排序的结果，这个结果相当于一个临时结果，所以我们不需要把它输出出来而可以通过linux下面的管道“|” (<- 这边的这个竖线就叫管道，pipeline)直接传递给下一个命令行，这边就是bcftools (后面我们会看到，还有其他的软件可以用来处理mpileup得到的结果)，所以bcftools最后一个参数是一个“-” ，表示标准输入(STDIN)，就是通过管道给它的东西(所以管道其实是个很形象的东西，他不需要你把中间的结果先保存下来再读取进去，而是直接连接两个命令行，把一个命令行的输出传递到另一个命令行的输入，当然如果你想的话，可以用管道一个连一个的连恩多个……)

 

关于mpileup的具体内容可以参见

http://samtools.sourceforge.net/mpileup.shtml

 

Samtools的优点是速度比较快，但是Call出来的Indel数量相比其他的好像稍微少一点。准确度方面不太好说，Samtools的主要开发者Heng Li先生(同样也是BWA、MAQ等等等等软件的主要开发者，超牛的说，瞬间民族自豪感max)自己也在给别人的回复中提到(我在BioStar的某个post里看到的= =)，对于不同的数据，其实每个软件的表现都不一致，有的软件适合某些数据，而有的软件就适合另外的数据，所以最重要的还是靠大家自己多发掘，找到适合自己数据的方法。不过Samtools停止开发已久，相比之下还是目前一直在开发的如GATK等等软件可能比较有优势一点(后来Heng Li先生也已经到broadinstitute去工作了OTL)。

 

 

(2) GATK UnifiedGenotyper

http://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_genotyper_UnifiedGenotyper.html

GATK之前的帖子已经讨论了很多了，call Indel和SNP的区别也不大，for example,

 

java -jar GenomeAnalysisTKLite.jar 
   -R ref.fasta 
   -T UnifiedGenotyper 
   -I sample.bam 
   -o sample.gatk.vcf 
   -nt 4 
   -stand_call_conf 50.0 
   -stand_emit_conf 0 
   -glm INDEL 
   -rf BadCigar

 

这边主要是把-glm这个参数设成了INDEL，所以输出的结果当中只有INDEL。

 

另外现在转用Lite版本了，终于可以不再受没有gatk_key带来的困扰了(- -||

 

-nt这个参数设的是线程数(我居然一直都不知道UnifiedGenotyper已经支持多线程了= =，现在可以设这个参数我们就再也不怕他跑得太慢啦……)，除了这个参数还有-nct也可以控制线程数，目测区别比较微妙，data threads是个什么概念楼主也比较费解，-nt是分配多少个data threads，而-nct是每个data threads分配多少个CPU，大家看机器资源试着设好了，我一般都用-nt，耗内存高一点不过应该貌似快一点……

 

然后根据个人经历2.2版本Call出来的Indel会莫名其妙的少掉很多，很多2.1 Call的出来的我用2.2试了恩个参数也Call不出来(肯定是我打开软件的方式不对= =)，而且2.2的AD值貌似有bug，明显很多不对，不过-maxAltAlleles默认值已经升到了6而且升值不减速确实很imba(就是这么多alt偶尔会觉得没什么意义的就是……) 不过如果真是bug就有望被修复，坐等GATK继续越做越好。

 

最后再提一下-rf这个参数，全称是--read_filter，就是用来筛选输入的bam文件中的reads的，因为GATK会检查bam文件里面有个叫Cigar值的东西，有时候有的mapping软件生成的bam文件当中有一些不符合它的标准，在用GATK处理时就可能会包Malformed read一类的错，所以可以通过-rf BadCigar这个参数来剔除掉这些不规范的reads，这样GATK就能正常运行了，上次有同学就碰到这样的问题，我后来才想起来加上这个参数应该大部分相关问题都能解决(如果加上了还不能解决的话那就可能是版本的bug了，GATK的论坛上貌似就有人碰到过这种情况，多换几个版本试试吧……)。

 

 

 

(3) Shore

http://www.1001genomes.org/software/shore.html

Ossowski, S. et al. Sequencing of natural strains of Arabidopsis thaliana with short reads. Genome Research 18, 2024–2033 (2008).

 

Shore是由做拟南芥的一群人开发的集mapping以及数据分析为一体的一套流程及软件。

引文点这个http://genome.cshlp.org/content/18/12/2024.full，国外超级大牛Detlef Weigel 注目，后来的1001基因组就基本上都是用的这一套(当然也还是他们一群人)。

相比与其他的软件，Shore使用起来流程要稍复杂些，因为他有自己的数据结构，所以先要进行数据的转换，然后对于paired-end他还有一个correct的步骤，默认只Call 1~3bp的Indel。

Shore其实是我最开始做mapping的时候用的一个软件(因为当初我也主要是做A.thaliana的= =)，但后来觉得真心不太好用，虽然貌似(有文献为证)还是很不错的一个软件(但目测就他们自己用，文章基本都是他们发的)，所以这边我就偷个懒直接跳过了，大家有兴趣的可以自行研究，功能真的其实还是很强大的……

 

 

(4) VarScan

http://varscan.sourceforge.net/

Koboldt, D. C. et al. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics 25, 2283–2285 (2009).

Koboldt, D. C. et al. VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).

 

引用率貌似还可以的一个variants检测软件，用来Call Indel自然不在话下。

前面我们说到samtools里面的mpileup，他生成的结果可以给bcftools用来Call Indel，当然也可以用其他软件来处理，这边的VarScan就也是通过mpileup的结果来Call Indel的，具体用法例如，

 

samtools mpileup -f ref.fasta sample.bam | 
   java -jar VarScan.v2.3.3.jar mpileup2indel 
   --output-vcf 1 
   --vcf-sample-list sample_names.list 
   > sample.varscan.vcf

 

VarScan是一个java程序，这--output-vcf 1表示输出结果格式为vcf格式的，否则就是软件本身的格式，然后--vcf-sample-list这个参数是可以不用加的，但是生成的vcf文件中sample名是按1、2、3…这样重新命名的，所以可以用这个参数给一个sample名称的列表，对应你给的bam文件中的sample名，这样vcf文件中就有对应的sample名称了。其他一些参数一般用默认的就ok，注意这边用的是mpileup2indel，对应samtools的mpileup，以前samtools里面还是pileup的时候他就对应pileup2indel (从我懂事起感觉pileup就被淘汰了，都没用到过，世事变迁啊……)，然后软件使用的时候可能会提示mpileup生成的结果里面有很多不能解析的，无视应该就可以了……

总体上来说使用方法还是很简单的，就相当于是把bcftools替换成了VarScan，相信大家很容易就能上手。

 

下面再介绍几个专门用于Call Indel的，软件名都是带Indel的，一看就很专业;-)

 

(5) Dindel

http://www.sanger.ac.uk/resources/software/dindel/

Albers, C. A. et al. Dindel: Accurate indel calls from short-read data. Genome Res. 21, 961–973 (2011).

 

又是sanger中心的东西，网上好像还能搜到另一个链接但貌似已经失效了(还是说被墙了= =)。

关于Dindel为什么要叫Dindel，个人猜测可能是detect Indel或者是discover Indel的简称，但不管是否正确，这个不是我们需要关心的问题…………

Dindel的输入文件也是bam文件，然后整个流程分四个步骤(好吧，又一个比较复杂的= =)，每一步的代码大致如下

 

## Stage 1 先把bam文件里面所有的Indel都提出来，在做这一步的同时软件
## 还会自动检测paired-end reads的insert size

dindel --ref ref.fasta 
   --analysis getCIGARindels 
   --bamFile sample.bam 
   --outputFile sample.dindel_output

第一步运行结束后会生成两个文件，我们这边就是

sample.dindel_output.variants.txt，

包含所有的候选indels，另外一个是

sample.dindel_output.libraries.txt，

里面是insert size的一个分布，下一步要用到的主要是第一个文件。

## Stage 2 把上面提出来的那些Indels划分成一个一个的windows(大小

## 为120bp左右)，这一步是用它提供的makeWindows.py这个python脚本

## 来执行的，输入文件就是第一步生成的sample.dindel_output.variants.txt

makeWindows.py --inputVarFile sample.dindel_output.variants.txt 
   --numWindowsPerFile 100000 
   --windowFilePrefix sample.dindel_output.windows

这一步根据你设置每个文件包含多少个windows即--numWindowsPerFile这个参数的大小，然后会生成几个前缀是sample.dindel_output.windows的文件，如
sample.dindel_output.windows.1.txt sample.dindel_output.windows.2.txt …..

## stage 3 这一步主要就是对每个windows里面的reads进行一个重排序，

## 这边要用到原始的bam文件，还要用到第一步生成的

## sample.dindel_output.libraries.txt这个文件，以及第二步生成的

## windows文件如sample.dindel_output.windows.1.txt


dindel --ref ref.fasta 
   --analysis indels 
   --doDiploid 
   --bamFile sample.bam 
   --varFile sample.dindel_output.windows.1.txt 
   --libFile sample.dindel_output.libraries.txt 
   --outputFile sample.dindel_output.stage2

然 后依然会生成好几个前缀为

sample.dindel_output.stage2

的文件，后缀名为*.glt.txt (这个软件不仅步骤多，生成的文件也好多啊- -)。这边大家了解一点的可能就已经发现，这边几步其实就跟以前我们讲的GATK做realign是一样的，这个告诉我们，在Indel周围做 realignment的步骤对于SNP和Indel准确性都是很重要的……

## Stage4 第四步就是生成最终结果的步骤

mergeOutputDiploid.py --ref ref.fasta 
   --inputFiles sample.dindel_output.stage2.list.txt 
   --outputFile sample.dindel_output.vcf

这边的输入文件sample.dindel_output.stage2.list.txt里面包含的是之前生成的所有glt.txt文件即sample.dindel_output.stage2.*.glt.txt文件的一个列表(不在当前路径下的话，文件名前面要带上相对路径，大家怕麻 烦的话就直接全用绝对路径好了)，然后会生成一个vcf的结果。

 

 

整个过程其实还算不太复杂，不过这边有一点要说的是Dindel的速度貌似有点不敢恭维，限速的主要是第三步，慢的当时我等了一天一夜才做完两三个samples实在看不下去然后毅然决然的把它掐了所以第四步我其实根本没跑过这种事我会告诉你？………… (一定是我的Indel比较难检测= +)，具体的还看大家自己的使用情况了。

 

 

 

(6) Pindel

https://trac.nbic.nl/pindel/

Ye, K., Schulz, M. H., Long, Q., Apweiler, R. & Ning, Z. Pindel: a pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads. Bioinformatics 25, 2865–2871 (2009).

 

跟其他Call Indel的软件不大一样，Pindel用的是一个叫pattern growth的算法来检测Indel以及其他的结构变异(所以才叫P-Indel的吧)，具体算法见上面的引文。引用次数还算可以，说明这种算法还是有一定优势的。

 

Pindel可以有几种输入文件方式，个人一般倾向于其中一种，而软件本身比较推荐的应该也是这种，具体流程代码其实很简单：

 

第一步就是Call Indel以及其他一些结构变异像倒位、串联重复等等

 

pindel -f ref.fasta -i sample.pindel.config -c ALL -T 2 -o sample

 

其实这一步我们就已经可以得到所有的结果，sample.pindel.config这个文件是一个配置文件，所有的bam文件以及insert size的信息就存放在这个文件里面，然后软件通过读取这个文件来作为它的输入，这个文件的内容格式如下

 

sample.bam  250 sample

 

第一列就是bam文件的文件名(需要时要带上路径)，第二列是insert size

 

* 这边补充讲一下insert size，所谓insert size就是你建库的时候序列打断的长度，就是你测的paired-end的序列就是从这样一条序列上面的两端，

 

Example:

|-----75----|------------------------100-----------------|-----75-----|

 

这边整个250bp其实就是你打断以后得到的序列的长度，然后两端各测了75bp，这个250bp就是我们讲的insert size，这个长度可以询问测序公司(没问过，不知道他们会不会告诉你。。。)，也可以通过软件来统计，像上面的Dindel的第一步就会先统计推测这个长度，然后picard里面也有CollectInsertSizeMetrics这个工具。

 
因为这个长度其实只是一个范围(一般比较窄)，打断长度基本都在这个范围里面，例如下图所示

 

 

所以这边只要设一个大概值就行，不用很精确；最后一列是设一个标签，因为这边可以设多个bam文件，这边的标签就会代替文件名出现在最终的结果中来区分reads的不同来源。列与列之间用制表符或者空格分开。

-c参数可以用来设定区域范围，-c ALL就表示整个基因组，-T是线程数，然后有个-w参数可以用来控制内存的使用量，大内存可以无视。

最后-o这个参数设的还是一个前缀，然后默认情况下会输出所有的插入缺失或者结构变异类型，分别生成以下后缀名结尾的文件：

 

D = deletion 缺失序列

SI = short insertion   短的插入序列

INV = inversion 转位

TD = tandem duplication      串联重复

LI = large insertion  长的插入序列，这个文件的格式跟其他文件的很不相同

BP = unassigned breakpoints       没有分到上面任意一种类型剩下来的断点

 

 

接下来我们可以通过里面提供的pindel2vcf这个程序把这些文件转换成我们常用的vcf文件，

方便下游处理，

pindel2vcf -r ref.fasta -R REF_NAME -d 20121208 -P sample -v sample.indel.vcf

 

 

这边-R需要为参考序列的设定一个名称，-d还要设定日期(就是这个参考序列生成的日期)，当然随便设应该也没什么问题，主要还是为了规范化，-P这个参数的设定值就是前面生成结果的那个前缀名例如这边是sample(当然所有的sample_*结果都得在这边)，最后-v是生成的vcf文件的文件名。

 

Pindel用起来还是比较简单的，而且速度还是比较快的，不过生成的vcf文件不是那么的标准，用GATK这种软件处理的时候可能不太方便，可以加上-G这个参数来让它尽可能符合GATK输入文件的要求。

 

 

(7) SOAPindel

http://soap.genomics.org.cn/soapindel.html

http://sourceforge.net/projects/soapindel/

Li, S. et al. SOAPindel: Efficient identification of indels from short paired reads. Genome Res. (2012).doi:10.1101/gr.132480.111

 

SOAP(Short Oligonucleotide Analysis Package)是由华大开发的一些列二代测序结果分析软件，当年我用SOAP的时候SOAPindel还处在自己人都不敢用的开发初期阶段，后来貌似SOAP团队的一些开发者都已经跳槽离开华大了(手头上还有他们公司给的传单……)，SOAP的网站也还一直挂着Under Construction For Full Updates却依旧保持这个风格……(刚发现现在有sourceforge页面了，大家也可以多留意关注)。

再后来，SOAPindel就出来了，其实出来也已经挺久的了，不过我还是看到Genome Research的文章的时候才想起来原来还有这个，BGI的软件貌似也是内部发的文章用的多(我口胡)…… 命令行其实也很简单，这边暂时就不多介绍了。等以后有机会用到再补上吧。