Pindel检测“复杂INDEL”的一次测试

2019-06-10 本文已影响2人徐广惠_6f76

写在前面

最近在实际工作中，遇到了一些“复杂InDEL”容易漏检的问题；这是由于在肿瘤组织中，一些插入和缺失可能同时发生于基因组的相同或相近的位置，造成了之前使用的SNV检测软件存在漏检风险。为了解决这个问题，这几天测试了Pindel软件在检测复杂INDEL方面的表现。

Pindel的使用

使用conda安装Pindel后，直接命令行运行pindel就可以查看软件的操作文档。其中基本的参数如下，如果要检测全部染色体，-c参数可以省略。

 pindel -f <reference.fa> -p <pindel_input>
           [and/or -i bam_configuration_file]
           -c <chromosome_name> -o <prefix_for_output_file>

其中第二个参数可以使用-p或-i，实话说我没有搞清楚-p参数需要输入的pindel_input文件到底是什么，所以选择了-i，bam_configuration_file是一个配置文件，所有的bam文件以及insert size的信息就存放在这个文件里面，然后软件通过读取这个文件来作为它的输入，这个文件的内容格式如下:

# bam_configuration_file用tab或空格分隔；第二列即为insert size，大概即可；可输入多个bam
sample.bam  150   sample

除了基本参数，pindel还有其他可选参数，其中我将-M参数由默认的1改为2，即变异支持数最少为2个，参数解释具体如下：

-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

但是这种结果文件的格式并不太易于阅读，也不方便后续处理，因为它看起来是这样的：

      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，可以将上述的结果文件转化为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 参数可以将所有结果文件转为VCF，-p参数可以指定其中的一个结果文件，上面示例就是只把Deletion的结果文件转为VCF。由于此次测试我只关注INDEL，所以只转出了两个文件:

-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数据分析用时不到半个小时。

第一个INDEL

Pindel的结果记录：

7    55242463        .       AAGG    A       .       PASS    END=55242466;HOMLEN=1;HOMSEQ=A;SVLEN=-3;SVTYPE=DEL     GT:AD   0/1:1259,876

第二个INDEL