寫在前面
最近在實際工作中，遇到了一些“復雜InDEL”容易漏檢的問題；這是由于在腫瘤組織中，一些插入和缺失可能同時發(fā)生于基因組的相同或相近的位置，造成了之前使用的SNV檢測軟件存在漏檢風險。為了解決這個問題，這幾天測試了Pindel軟件在檢測復雜INDEL方面的表現(xiàn)。
Pindel的使用
使用conda安裝Pindel后，直接命令行運行pindel就可以查看軟件的操作文檔。其中基本的參數(shù)如下，如果要檢測全部染色體，-c參數(shù)可以省略。
pindel -f <reference.fa> -p <pindel_input>
[and/or -i bam_configuration_file]
-c <chromosome_name> -o <prefix_for_output_file>
其中第二個參數(shù)可以使用-p或-i，實話說我沒有搞清楚-p參數(shù)需要輸入的pindel_input文件到底是什么，所以選擇了-i，bam_configuration_file是一個配置文件，所有的bam文件以及insert size的信息就存放在這個文件里面，然后軟件通過讀取這個文件來作為它的輸入，這個文件的內(nèi)容格式如下:
# bam_configuration_file用tab或空格分隔；第二列即為insert size，大概即可；可輸入多個bam
sample.bam 150 sample
除了基本參數(shù)，pindel還有其他可選參數(shù)，其中我將-M參數(shù)由默認的1改為2，即變異支持數(shù)最少為2個，參數(shù)解釋具體如下：
-M/--minimum_support_for_event
Pindel only calls events which have this number or more supporting
reads (default 1)
結果文件
Pindel將不同的SV輸出在不同的結果文件中，并在文件名中標注變異類型:
-rw-r--r-- 1 xu informatics 0 Jun 6 14:32 sample_BP # BP = unassigned breakpoints
-rw-r--r-- 1 xu informatics 0 Jun 6 14:32 sample_CloseEndMapped
-rw-r--r-- 1 xu informatics 259655 Jun 6 14:34 sample_D # D = deletion
-rw-r--r-- 1 xu informatics 0 Jun 6 14:35 sample_INT_final
-rw-r--r-- 1 xu informatics 8037 Jun 6 14:33 sample_INV # INV = inversion
-rw-r--r-- 1 xu informatics 0 Jun 6 14:32 sample_LI # LI = large insertion
-rw-r--r-- 1 xu informatics 0 Jun 6 14:32 sample_RP
-rw-r--r-- 1 xu informatics 242141 Jun 6 14:34 sample_SI # SI = short insertion
但是這種結果文件的格式并不太易于閱讀，也不方便后續(xù)處理，因為它看起來是這樣的：
1 ####################################################################################################
2 0 D 1 NT 0 "" ChrID 1 BP 20915589 20915591 BP_range 20915589 20915593 Supports 11 11 + 6 6 - 5 5 S1
3 GAGTCCAAACCATGGGAGGCTCCTCTCCTAGACCCTGCATCCTGAAAGCTGCGTACCTGAGAGCCTGCGGTCTGGCTGCAGGGACACACCCAAGGGGAGGAGCTGCAATCGTGTCTGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTTcCCGCTGCTCTG
4 CCTGTTT CCGCTGCTCTG
5 CCTGTTT CCGCTGCTCTG
6 CCTGTTT CCGCTGCTCTG
7 TCCTGTTT CCGCTGCTCTG
8 TCCTGTTT CCGCTGCTCTG
9 TCCTGTTT CCGCTGCTCTN
10 TGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTT CCGCTGCTCTG
11 CTGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTT CCGCTGCTCTG
12 TGCAATCGTGTCTGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTT CCGCTGCTCTN
13 AGGAGCTGCAATCGTGTCTGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTT CCGCTGCTCTG
14 CTGAGAGCCTGCGGTCTGGCTGCAGGGACACACCCAAGGGGAGGAGCTGCAATCGTGTCTGGGGCCCCAGCCCAGGCTGGCCGGAGCTCCTGTTT CCGCTGCTCTG
15 ####################################################################################################
16 1 D 6 NT 0 "" ChrID 1 BP 26235085 26235092 BP_range 26235085 26235136 Supports 6 6 + 6 6 - 0 0 S1
17 TGGGTTTGTTGGGAGATGCCTGGTTCTGCCAGTCTGTGACAATGTTCCAAGCTCCTCACAGCTGCTTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTTttattaTTATTA
18 ATGAGCCTTTGAATAANAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
19 CCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
20 TAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
21 GGCTGGGTTNAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
22 TTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
23 GCTCCTCACAGCTGCTTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATTCCTGTCCTT TTATTA
24 ####################################################################################################
25 2 D 19 NT 13 "GATTCCTGTCCTT" ChrID 1 BP 26235075 26235095 BP_range 26235075 26235095 Supports 2 2 + 2 2 - 0
26 ACTTGGCTTCTGGGTTTGTTGGGAGATGCCTGGTTCTGCCAGTCTGTGACAATGTTCCAAGCTCCTCACAGCTGCTTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGAT
27 TTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATGATTCCTGTCCT
28 TGCTTGAGGACTGAGAGGGCTGGGTTAAAGTTTCCCCTAGAATGAGCCTTTGAATAAAAAGGTGCTTTTGAGGTGGGATGATTCCTGTCCT
29 #########
因此Pindel很貼心的提供了一個腳本——pindel2vcf，可以將上述的結果文件轉(zhuǎn)化為vcf格式。pindel2vcf使用的方法如下：
pindel2vcf -p sample3chr20_D -r human_g1k_v36.fasta -R 1000GenomesPilot-NCBI36
-d 20101123 -v sample3chr20_D.vcf
or (with -P): pindel2vcf -P sample3chr20 -r human_g1k_v36.fasta -R 1000GenomesPilot-NCBI36 -d 20101123 -v sample3chr20_all.vcf
使用-P 參數(shù)可以將所有結果文件轉(zhuǎn)為VCF，-p參數(shù)可以指定其中的一個結果文件，上面示例就是只把Deletion的結果文件轉(zhuǎn)為VCF。由于此次測試我只關注INDEL，所以只轉(zhuǎn)出了兩個文件:
-rw-r--r-- 1 xu informatics 4352 Jun 6 14:40 sample_del.vcf
-rw-r--r-- 1 xu informatics 9079 Jun 6 14:43 sample_insertion.vcf
測試結果
作為測試的兩個復雜INDEL最終都被Pindel成功檢出了，100M左右的bam數(shù)據(jù)分析用時不到半個小時。
第一個INDEL
Pindel的結果記錄：
7 55242463 . AAGG A . PASS END=55242466;HOMLEN=1;HOMSEQ=A;SVLEN=-3;SVTYPE=DEL GT:AD 0/1:1259,876
第二個INDEL
Pindel的結果記錄：
7 55242467 . AATTAAGAGAAG AGC . PASS END=55242478;HOMLEN=0;SVLEN=-11;SVTYPE=RPL;NTLEN=2 GT:AD 0/1:327,137
