selscan dropping variants?

[mmcqui]

I am trying to understand how selscan excludes variants, as a region I am interested in appears to be missing some important variants in my output. My command is:

selscan --ihs --vcf inFile.vcf.gz --map mapFile.map --out population.iHS --threads 25

My input .vcf file contains 4,702,160 variants (many of which are rare). The output .log file tells me it dropped 4,298,509 variants due to low frequency (maf < 0.05), which is correct based on my external checks. The .log file also tells me it dropped an additional 10,622 variants due to large gaps, reaching the chromosome edge, etc. This should leave us with 393,029 variants in the iHs .out file. However, the .out file contains slightly fewer, 392,297 variants. Is there somewhere else that selscan is dropping variants, that isn't obvious from the .log file? Perhaps due to the EHH cutoff value, which is 0.05 by default?

Thanks for any information you can provide.

[mmcqui]

After looking at some 'closed' issues on this site, I found someone with a similar issue. I realized that many of my variants in my region of interest were assigned the same genetic position in the genetic map file, which I'm guessing is why they were dropped?. In case others have this issue working with human data, the widely available genetic maps based on Hapmap contain large regions where all variants have the same genetic position. If you're doing any kind of interpolation for your variants to include in selscan, it appears the program will drop these variants from your output.

[szpiech]

Hi Michael, Usually I linearly interpolate genetic map values for sites in between positions for which I have map information, there's a simple program in my github that can do this (just mind the --max-gap flag, the default is quite small). Usually if the precision is high enough it doesn't create too many issues with sites being assigned exactly the same position.

…

-Zachary

On Wed, Aug 21, 2024 at 12:44 PM Michael McQuillan ***@***.***> wrote: After looking at some 'closed' issues on this site, I found someone with a similar issue. I realized that many of my variants in my region of interest were assigned the same genetic position in the genetic map file, which I'm guessing is why they were dropped?. In case others have this issue working with human data, the widely available genetic maps based on Hapmap contain large regions where all variants have the same genetic position. If you're doing any kind of interpolation for your variants to include in selscan, it appears the program will drop these variants from your output. — Reply to this email directly, view it on GitHub <#124 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABAKRQXUCRZVW4F7NN2ZZNLZSS7XZAVCNFSM6AAAAABMX4E5B2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBSGUZTEMJUGY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[mmcqui]

Thanks for the response! Yes, I currently use linear interpolation for sites in between. The issue I think is that I am using publicly available genetic maps for my 'known' positions, and these maps have many sites with identical genetic positions, causing my linear interpolation to also assign any variant within these regions the same genetic map position. Seems like I'll need to find new maps.

[szpiech]

Ahh, ok. I think a simple solution could be to drop sites from your reference map such that only sites with unique genetic positions are included, then do linear interpolation. Zachary Le mer. 21 août 2024 à 13:30, Michael McQuillan ***@***.***> a écrit :

…

Thanks for the response! Yes, I currently use linear interpolation for sites in between. The issue I think is that I am using publicly available genetic maps for my 'known' positions, and these maps have many sites with identical genetic positions, causing my linear interpolation to also assign any variant within these regions the same genetic map position. Seems like I'll need to find new maps. — Reply to this email directly, view it on GitHub <#124 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABAKRQTVCNAKWXIY32MJVL3ZSTFCTAVCNFSM6AAAAABMX4E5B2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBSGYYTAOJUHE> . You are receiving this because you commented.Message ID: ***@***.***>

[mmcqui]

Ah, great idea! Thanks so much!