search

GeneDx / scramble

Other
36 stars 12 forks source link

Help with call filtering #7

Closed hackdna closed 4 years ago

hackdna commented 4 years ago

We ran SCRAMble on NA12878 dataset and compared the output with NA12878 validation data.

Could you help us explain the remaining differences in calls? Could you recommend any additional filters to use for calls made by SCRAMble?

rebeccaito commented 4 years ago

Thanks for your question. If you want to increase the specificity of Scramble, you can use the -n option to increase the minimum number of clipped reads and the --pct-align option to increase the threshold for percent length of clipped read that aligns to MEI sequence. You can also perform post hoc filtering on the Scramble output file. For example, since most Alu MEIs are full length or mostly full length, you can require that Alus have Start_in_MEI that is < 20. This will remove MEI calls that might be false positives due to other kinds of structural variation that involve Alu sequence at the breakpoints. You can also filter on Alignment_Percent_Identity as lower values are more likely to be false positives. If you do filter on Alignment_Percent_Identity, you may want to set different thresholds for Alu, L1, and SVA since SVAs can have lower Alignment_Percent_Identity and still be real calls.

I will note that the ideal Scramble use-case is for targeted sequencing data. Was this run on WGS? Also, the false negative calls identified by MELT, might be due to MEIs with 5' transductions which we would not expect Scramble to detect.

Happy to discuss further!

hackdna commented 4 years ago

Thank you for a quick and informative reply. I have applied Alu filtering using Start_in_MEI < 20 and it has reduced both the number of false positives and real calls: index

Positives and negatives by class: pos_neg_by_class

Alignment_Percent_Identity distributions per MEI class for true and false positives: index index

Also, I was wondering why is targeted sequencing an ideal use case? Is it just depth?

rebeccaito commented 4 years ago

Targeted sequencing is less likely to overlap with repetitive sequence in the genome and is less likely to produce false positive clipped clusters. Also, for performance, Scramble was intended for targeted NGS and may have memory issues with WGS.

Do you have a specific sensitivity/specificity goal for your use-case when running Scramble? Have you done a deep dive to characterize the false positive/negative calls? For example, MEIs with 5' transductions are not expected to picked up by Scramble.

rebeccaito commented 4 years ago

Also, I would recommend using the PALMER MEI calls from the GIAB data set as a gold standard.

hackdna commented 4 years ago

Thank you very much for your suggestions. I just wanted to share what we got when comparing with PALMER calls (obtained from https://academic.oup.com/nar/article/48/3/1146/5680708): scramble-palmer

rebeccaito commented 4 years ago

You might be using different bams than we did for the benchmarking in the Genetics in Medicine paper (https://www.nature.com/articles/s41436-020-0749-x, supporting info page 4). We had a sensitivity of 85% at default settings.

We downloaded the GIAB sample from ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG002_NA24385_son/NIST_HiSeq_HG002_Homogeneity-10953946/HG002Run01-11419412/HG002run1_S1.bam, and ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG002_NA24385_son/NIST_HiSeq_HG002_Homogeneity-10953946/HG002Run02-11611685/HG002-Run2_S1.bam).

Do you have a specific question I can help with? Can you provide a reproducible example of the issue you are having?

hackdna commented 4 years ago

I've used BAMs produced by this workflow: https://bcbio-nextgen.readthedocs.io/en/latest/contents/germline_variants.html#workflow4-whole-genome-trio-50x-hg38 We simply wanted to see how well Scramble performs before adding it to bcbio. The remaining step is to make Scramble available via Bioconda (#6). Thanks again for all your help.