How to calculate frequencies of SNPs with respect to a particular allele by plink?
1
0
Entering edit mode
3.1 years ago
jamespower ▴ 100

Is there a way to force Plink to compute frequency with respect to a specific allele (NOT minnor allele frequency), for example using "--a2-allele" with a list of reference alleles?

plink -bfile --freq --extract list.snps --a2-allele list.reference.allele

Thank you

frequency snp plink • 3.0k views
ADD COMMENT
1
Entering edit mode

Hi James, if you receive no answer here, then I would suggest that you post this question in the PLINK2 user group: https://groups.google.com/g/plink2-users

ADD REPLY
2
Entering edit mode
3.1 years ago

Your proposed approach should work.

ADD COMMENT
0
Entering edit mode

Hi there, the approach I have computes the MAF not the frequency, while I would like the frequency based on the specified allele.

ADD REPLY
0
Entering edit mode

I just tried this, and some of the reported frequencies were above 0.5. So this is definitely not just reporting MAFs. Please elaborate, with a full .log file and an example of an unexpected result.

ADD REPLY
0
Entering edit mode

Hi, When I try, the frequencies are always max 0.5 (reproducible example and log below). You can clearly see it ignores the option because when flipping the command to "--a1-allele list.reference.allele" instead of "--a2-allele list.reference.allele", nothing changes...

ex. list.snps

rs139938544
rs10841218
rs12817774
rs11044517
rs57357758
rs79072570
rs2115857
rs11044518
rs148808928
rs2961367
rs66703594
rs16915433
rs77446069
rs150097857
rs149238240
rs117083095
rs140136853
rs3752823
rs8601
rs843869

ex. list.reference.allele

rs139938544     t
rs10841218      a
rs12817774      t
rs11044517      t
rs57357758      a
rs79072570      a
rs2115857       t
rs11044518      t
rs148808928     t
rs2961367       t
rs66703594      a
rs16915433      t
rs77446069      t
rs150097857     a
rs149238240     t
rs117083095     a
rs140136853     a
rs3752823       a
rs8601  t
rs843869        t

This is the command and log file:

plink -bfile EUR/chr12 --freq --extract list.snps --a2-allele list.reference.allele --out snps_freq

Log file reported here:

PLINK v1.90b6.10 64-bit (17 Jun 2019)
Options in effect:
  --a2-allele list.reference.allele
  --bfile EUR/chr12
  --extract list.snps
  --freq
  --out snps_freq

Start time: Thu Mar 25 11:38:35 2021

Random number seed: 1616668715
64164 MB RAM detected; reserving 32082 MB for main workspace.
3848775 variants loaded from .bim file.
503 people (0 males, 0 females, 503 ambiguous) loaded from .fam.
Ambiguous sex IDs written to
snps_freq.nosex
.
--extract: 7354 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 503 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Warning: Impossible A2 allele assignment for variant rs139938544.
Warning: Impossible A2 allele assignment for variant rs10841218.
Warning: Impossible A2 allele assignment for variant rs12817774.
Warning: Impossible A2 allele assignment for variant rs11044517.
Warning: Impossible A2 allele assignment for variant rs57357758.
....
--a2-allele: 7354 assignments made.
--freq: Allele frequencies (founders only) written to
snps_freq.frq
ADD REPLY
0
Entering edit mode

Your list.reference.allele has lowercase alleles, and the "Impossible A2 allele assignment" warnings make me believe that your main dataset has uppercase alleles.

I will fix the --a2-allele log message in the next plink 1.9 build to report that no allele assignments were actually made here because there were mismatches for all the allele codes.

ADD REPLY
0
Entering edit mode

Hi chrchang, ok! thank you! I confirm this works when I properly match the uppercase. So to be clear the output for .frq is the MAF, and not the frequency, because the --a2 option is ignored since the data in list.reference.allele and the data in the genotypes cannot be matched? In case this is the output created, it would be great to not create an output with the MAF, but a warning instead, to be able to catch these cases. Thank you much for your help.

ADD REPLY

Login before adding your answer.

Traffic: 1823 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6