Change logs:

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

Indel Calling相比于SNP Calling的難度要大一些，因?yàn)橛捎谶@種插入-缺失的存在,，本身就很容易干擾排序,，這種干擾會導(dǎo)致Indel周圍出現(xiàn)很多假陽性的SNP，而且會影響Indel本身的準(zhǔn)確性,。理論上來說,，檢測Indel的最好方式就是做de novo assembly，然后比較de novo得到的基因組與原來的基因組,，不過實(shí)際上de novo assembly的難度更大OTL

Paired-end測序?yàn)閷ふ逸^長片段的Indel提供了非常有用的信息,，但是如何準(zhǔn)確的利用這些信息也是目前這一塊的一大難點(diǎn)。

下面就簡單的介紹一些現(xiàn)有的Call Indel的軟件及使用流程,，注意一下這邊介紹的幾個基本上都是針對Illumina的paired-end數(shù)據(jù)的,。

(1) Samtools mpileup

http://samtools./

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

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

這邊有個-N的參數(shù)表示跳過reference里面堿基是N的位點(diǎn),。

注意這邊其實(shí)用的是一個管道，因?yàn)?/span>mpileup直接產(chǎn)生的結(jié)果不是最后vcf文件里看到的結(jié)果,，而是一個排序的結(jié)果,，這個結(jié)果相當(dāng)于一個臨時(shí)結(jié)果，所以我們不需要把它輸出出來而可以通過linux下面的管道“|” (<- 這邊的這個豎線就叫管道,，pipeline)直接傳遞給下一個命令行,，這邊就是bcftools (后面我們會看到，還有其他的軟件可以用來處理mpileup得到的結(jié)果),，所以bcftools最后一個參數(shù)是一個“-” ,，表示標(biāo)準(zhǔn)輸入(STDIN)，就是通過管道給它的東西(所以管道其實(shí)是個很形象的東西,，他不需要你把中間的結(jié)果先保存下來再讀取進(jìn)去，而是直接連接兩個命令行,，把一個命令行的輸出傳遞到另一個命令行的輸入,，當(dāng)然如果你想的話，可以用管道一個連一個的連恩多個……)

關(guān)于mpileup的具體內(nèi)容可以參見

http://samtools./mpileup.shtml

Samtools的優(yōu)點(diǎn)是速度比較快,，但是Call出來的Indel數(shù)量相比其他的好像稍微少一點(diǎn),。準(zhǔn)確度方面不太好說，Samtools的主要開發(fā)者Heng Li先生(同樣也是BWA,、MAQ等等等等軟件的主要開發(fā)者,，超牛的說，瞬間民族自豪感max)自己也在給別人的回復(fù)中提到(我在BioStar的某個post里看到的= =),，對于不同的數(shù)據(jù),，其實(shí)每個軟件的表現(xiàn)都不一致,，有的軟件適合某些數(shù)據(jù)，而有的軟件就適合另外的數(shù)據(jù),，所以最重要的還是靠大家自己多發(fā)掘,，找到適合自己數(shù)據(jù)的方法。不過Samtools停止開發(fā)已久,，相比之下還是目前一直在開發(fā)的如GATK等等軟件可能比較有優(yōu)勢一點(diǎn)(后來Heng Li先生也已經(jīng)到broadinstitute去工作了OTL),。

(2) GATK UnifiedGenotyper

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

GATK之前的帖子已經(jīng)討論了很多了，call Indel和SNP的區(qū)別也不大,，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這個參數(shù)設(shè)成了INDEL,，所以輸出的結(jié)果當(dāng)中只有INDEL。

另外現(xiàn)在轉(zhuǎn)用Lite版本了,，終于可以不再受沒有gatk_key帶來的困擾了(- -||

-nt這個參數(shù)設(shè)的是線程數(shù)(我居然一直都不知道UnifiedGenotyper已經(jīng)支持多線程了= =,，現(xiàn)在可以設(shè)這個參數(shù)我們就再也不怕他跑得太慢啦……)，除了這個參數(shù)還有-nct也可以控制線程數(shù),，目測區(qū)別比較微妙,，data threads是個什么概念樓主也比較費(fèi)解，-nt是分配多少個data threads,，而-nct是每個data threads分配多少個CPU,，大家看機(jī)器資源試著設(shè)好了，我一般都用-nt,，耗內(nèi)存高一點(diǎn)不過應(yīng)該貌似快一點(diǎn)……

然后根據(jù)個人經(jīng)歷2.2版本Call出來的Indel會莫名其妙的少掉很多,，很多2.1 Call的出來的我用2.2試了恩個參數(shù)也Call不出來(肯定是我打開軟件的方式不對= =)，而且2.2的AD值貌似有bug,，明顯很多不對,，不過-maxAltAlleles默認(rèn)值已經(jīng)升到了6而且升值不減速確實(shí)很imba(就是這么多alt偶爾會覺得沒什么意義的就是……) 不過如果真是bug就有望被修復(fù)，坐等GATK繼續(xù)越做越好,。

最后再提一下-rf這個參數(shù),，全稱是--read_filter，就是用來篩選輸入的bam文件中的reads的,，因?yàn)?/span>GATK會檢查bam文件里面有個叫Cigar值的東西,，有時(shí)候有的mapping軟件生成的bam文件當(dāng)中有一些不符合它的標(biāo)準(zhǔn)，在用GATK處理時(shí)就可能會包Malformed read一類的錯,，所以可以通過-rf BadCigar這個參數(shù)來剔除掉這些不規(guī)范的reads,，這樣GATK就能正常運(yùn)行了，上次有同學(xué)就碰到這樣的問題,，我后來才想起來加上這個參數(shù)應(yīng)該大部分相關(guān)問題都能解決(如果加上了還不能解決的話那就可能是版本的bug了,，GATK的論壇上貌似就有人碰到過這種情況，多換幾個版本試試吧……),。

(3) Shore

http://www./software/shore.html

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

Shore是由做擬南芥的一群人開發(fā)的集mapping以及數(shù)據(jù)分析為一體的一套流程及軟件,。

引文點(diǎn)這個http://genome./content/18/12/2024.full,，國外超級大牛Detlef Weigel 注目，后來的1001基因組就基本上都是用的這一套(當(dāng)然也還是他們一群人),。

相比與其他的軟件,，Shore使用起來流程要稍復(fù)雜些，因?yàn)樗凶约旱臄?shù)據(jù)結(jié)構(gòu),，所以先要進(jìn)行數(shù)據(jù)的轉(zhuǎn)換,，然后對于paired-end他還有一個correct的步驟，默認(rèn)只Call 1~3bp的Indel,。

Shore其實(shí)是我最開始做mapping的時(shí)候用的一個軟件(因?yàn)楫?dāng)初我也主要是做A.thaliana的= =),，但后來覺得真心不太好用，雖然貌似(有文獻(xiàn)為證)還是很不錯的一個軟件(但目測就他們自己用,，文章基本都是他們發(fā)的),，所以這邊我就偷個懶直接跳過了，大家有興趣的可以自行研究,，功能真的其實(shí)還是很強(qiáng)大的……

(4) VarScan

http://varscan./

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,，他生成的結(jié)果可以給bcftools用來Call Indel,，當(dāng)然也可以用其他軟件來處理，這邊的VarScan就也是通過mpileup的結(jié)果來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表示輸出結(jié)果格式為vcf格式的,，否則就是軟件本身的格式,，然后--vcf-sample-list這個參數(shù)是可以不用加的，但是生成的vcf文件中sample名是按1,、2,、3…這樣重新命名的，所以可以用這個參數(shù)給一個sample名稱的列表,，對應(yīng)你給的bam文件中的sample名,，這樣vcf文件中就有對應(yīng)的sample名稱了。其他一些參數(shù)一般用默認(rèn)的就ok,，注意這邊用的是mpileup2indel，對應(yīng)samtools的mpileup,，以前samtools里面還是pileup的時(shí)候他就對應(yīng)pileup2indel (從我懂事起感覺pileup就被淘汰了,，都沒用到過，世事變遷啊……),，然后軟件使用的時(shí)候可能會提示mpileup生成的結(jié)果里面有很多不能解析的,，無視應(yīng)該就可以了……

總體上來說使用方法還是很簡單的,，就相當(dāng)于是把bcftools替換成了VarScan，相信大家很容易就能上手,。