加载中…
个人资料
邵先生
邵先生
  • 博客等级:
  • 博客积分:0
  • 博客访问:60,922
  • 关注人气:35
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

[转载]Indel Calling (Part2)

(2013-02-18 14:25:34)
标签:

转载

原文地址:Indel Calling (Part2)作者:Nowind

下面再介绍几个专门用于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 readsinsert 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 (这个软件不仅步骤多,生成的文件也好多啊- -)。这边大家了解一点的可能就已经发现,这边几步其实就跟以前我们讲的GATKrealign是一样的,这个告诉我们,在Indel周围做 realignment的步骤对于SNPIndel准确性都是很重要的……

## 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这个工具。

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

 [转载]Indel <wbr>Calling <wbr>(Part2)



 

所以这边只要设一个大概值就行,不用很精确;最后一列是设一个标签,因为这边可以设多个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的软件貌似也是内部发的文章用的多(我口胡)…… 命令行其实也很简单,这边暂时就不多介绍了。等以后有机会用到再补上吧。

 

 

 

 

 

后记,

1)      word里敲完复制到blog里面有的排版就乱了,让人有点淡淡的忧伤……

2)      我刚刚一点发布就只剩一个标题和前两行了是神马个情况!!!= =



你妹的居然是因为超了字符数限制,怎么从来没人告诉我还有字符数限制(泪流满面),png的图片又不行,这是要闹哪样……

0

  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有