Help with call filtering about scramble HOT 7 CLOSED

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!

from scramble.

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:

Positives and negatives by class:

Alignment_Percent_Identity distributions per MEI class for true and false positives:

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

from scramble.

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.

from scramble.

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

from scramble.

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):

from scramble.

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?

from scramble.

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.

from scramble.