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
(2) GATK UnifiedGenotyper
參考別人的一個流程：
http://blog.sina.com.cn/s/blog_681aaa550101bhdc.html
http://blog.sina.com.cn/s/blog_6721167201018fyw.html
http://blog.sina.com.cn/s/blog_6721167201018jik.html
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
內容來自於：https://i.cnblogs.com/EditPosts.aspx?opt=1