search

NCGG-MGC / IMSindel

IMSindel: An accurate intermediate-size indel detection tool incorporating de novo assembly and gapped global-local alignment with split read analysis
https://www.nature.com/articles/s41598-018-23978-z
MIT License
15 stars 0 forks source link

False negative homozygous duplication 77bps (or coordinates issue) #20

Open husamia opened 2 years ago

husamia commented 2 years ago

I have illumina paired reads 150bps that is targeted. The average coverage is pretty high. There is a known duplication of 77bps that is homozygous. It looks very good in IGV see snapshot.

image

you can see clipped reads on both ends. So I ran IMSindel and it detected several INS but non of them seems to be in the location of this known duplication. Also the depth seems to be off. Here is the commands:

docker command docker run -it --rm -v /mnt/e/:/home -w /home davelabhub/imsindel:latest bash

the command to run time imsindel --bam 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam --chr 1 --outd 21GN-187-0021/imsindel_out --indelsize 80 --reffa ref/human_g1k_v37.fasta --thread 10 --temp 21GN-187-0021/temp

log

`

Parameters: Avg. base quality: 20 Maping quality: 20 Read group: within 3bp paired B and F: within 5bp Support reads for making consensus sequence: 3 mimimum clipping fragment base: 5bp support clip length: 5bp bam: 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam chr: 1 outd: 21GN-187-0021/imsindel_out indelsize: 80 reffa: ref/human_g1k_v37.fasta glsearch: glsearch36 glsearch mat: /opt/IMSindel/lib/../data/mydna.mat mafft: mafft samtools: samtools temp: 21GN-187-0021/temp thread: 25 exclude-region

  1. collecting indel related reads... samtools view -F 1024 -f 2 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam 1

    backward_clips: 177649

    forward_clips: 177642

    non_clips: 18247

  2. collecting indel related reads...done
  3. collecting unmapped reads... samtools view -F 1024 -f 8 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam 1 mate_unmapped_read_names: 1063 samtools view -F 1024 -f 4 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam 1 Insert size Avg: 129.54610340127292 SD: 36.41371592197686

    unmapped reads: 460

  4. collecting unmapped reads...done
  5. considering support reads...

    backward clip with support reads: 2974

    forward clip with support reads: 2982

    non_clips with suport reads: 1298

  6. considering support reads...done
  7. making consensus seqs from support reads...

    backward clip with consensus: 2974 --> 2355

    forward clip with consensus: 2982 --> 2367

    shot indel with consensus: 1298 --> 154

  8. making consensus seqs from support reads...done
  9. making consensus seq from B and F.. making consensus seq for long deletion...done
  10. detection of indels...

    paired long indel candidates: 12

    unpaired long indel candidates: 464

    short indel candidates: 154

    all_pos: 10

    fin_indels: 10

    samtools view -F 1024 -f 2 21GN-187-0021/21GN-187-0021-MT_EB2_S13.bam 1

    all_pos_depth 10

  11. detection of indels...done

real 418m45.794s user 3040m11.668s sys 50m33.437s `

the output

indel_type call_type chr sttpos endpos indel_length INS,INS Hete 1 38,511,568 38,511,568 79,51 INS Hete 1 121,484,915 121,484,915 71 INS Hete 1 183,201,508 183,201,508 85 INS Hete 1 183,201,584 183,201,584 28 INS Hete 1 209,796,804 209,796,804 72 INS,INS Hete 1 209,797,351 209,797,351 37,93 INS Hete 1 215,822,056 215,822,056 78 INS Hete 1 215,822,181 215,822,181 37 INS Hete 1 216,373,025 216,373,025 41 INS Hete 1 216,373,037 216,373,037 70

the duplication is described as

LAMB3 (Laminin subunit beta 3) Variation Duplication (77 bps) in exon 10. This variation creates a frame shift starting at codon Asn345. The new reading frame ends in a STOP codon at position 77. This variant is known to ClinVar (August-2021): RCV000318696.6 (Pathogenic* - not provided), RCV000588290.1 (Pathogenic - Junctional epidermolysis bullosa gravis of Herlitz). This variant is known to gnomAD (2.1) <Exomes+Genomes>: ALL:0.0042% - AFR:0.040% - AMR:0.0056% HGVS Nomenclature (v15.11)

cDNA Level: NM_000228.2:c.958_1034dup gDNA Level: Chr1(GRCh37):g.209803180_209803256dup Protein Level: p.(Asn345Lysfs*77)

I am trying to run it with indel size 10,000 now. It would be nice to have VCF output so I can load it in IGV. The coordinates seems to be way off. I checked those locations and I don't see soft clipper reads so I am not sure what's going on.

holrock commented 2 years ago

Thank you for report. Since the IMSindel cannot detect all of intermediate indels, it could be false negative.