General Question About QTL Effect

[zjwinn]

Hi again,

This is more of a question rather than a bug or comment on a new feature. I am trying to pull the "true" effect of QTL out of the simulation parameters. I do this by going into the SP object and this allows me to pull something out of the slot called "addEff" which I am assuming are the additive effects of the QTL assigned for the trait. I am trying to compare this to GWAS estimate results. I am wondering, how is the effect calculated and what allele is considered the reference allele? I am assuming that the effect is calculated in reference to the major allele (e.g., 0) and the effect represents the additive effect of inheriting a copy of the alternate allele.

So, If I am comparing the estimated effect of an MTA identified by GWAS and the "true" additive effect of the QTL, then I need to make sure the effect is always calculated in reference to the reference allele (0) and expressed as inheriting the alternate allele (2).

Is this assumption correct?

Thanks,
Zach

[gregorgorjanc]

Allele 0 in AlphaSimR is assumed to the the ancestral allele if you have used runMacs(). If you have used quickHaplo() then the 0 and 1 have no real meaning. At any rate, which allele is 0 and 1 does not matter in terms of the additive effect - if you swap the 0 and 1 you just have to swap the sign of the additive effect;)

If you are simulating just an additive trait, then "addEff" is the thing you are after for comparing with GWAS, but note that "addEff" is a simulated effect (pulled from the assumed trait model), hence it's assumed known, and therefore does not need to be calculated. genParam() does return theoretical regression coefficients (alpha, allele substitution effects), which are akin to what you get from GWAS, but it really depends on what GWAS model you are fitting. Also, these true effects are calculated one at a time without incurring linkage&LD, so it's a theoretical quantity.

As to comparing between AlphaSimR and GWAS, yes you should be mindful how GWAS software you use will treat 0 and 1 (also be mindful that some software might consider 1 as major if it's more frequent). Do you really need to focus on the sign though? You could just look at the magnitude of the absolute effects too!?

@gaynorr another one to be moved to discussion forum;)

[zjwinn]

Oh yes, now that I see the discussion tab, this would definitely fit in that category. So, if I am pulling the "addEff" slot from "SP" for the specific trait, then that is the additive effect in either direction. I did start with runMacs() so I will go forth with the assumption that 0 is the ancestral allele and 1 is the alternate allele. It is true that I could look at the absolute value of the effect and compare the magnitudes. I think I will go forward with this assumption.

[gaynorr]

# Example GWAS analysis for K model using rrBLUP
library(AlphaSimR) 
library(rrBLUP) #For GWAS
library(qqman) #For Manhattan and Q-Q plots

# Global parameters
nQtlPerChr = 10
nSnpPerChr = 500 #must be >= nQtlPerChr
H2 = 0.8

#Simulate a mapping population
founderPop = runMacs(nInd=1000,
                     nChr=10,
                     segSites=nSnpPerChr,
                     inbred=TRUE)

# Setting simulation parameters
# Changing overlap for QTL and SNPs to true so that some SNPs will have an actual effect
SP = SimParam$
  new(founderPop)$
  restrSegSites(overlap=TRUE)$
  addTraitA(nQtlPerChr)$
  addSnpChip(nSnpPerChr)$
  setVarE(H2=H2)

pop = newPop(founderPop)

# Format genotypes for rrBLUP
rawGeno = pullSnpGeno(pop)
rawGeno = rawGeno-1
geno = as.data.frame(t(rawGeno))
geno = data.frame(snp=row.names(geno),
                  chr=rep(1:10,each=nSnpPerChr),
                  pos=rep(1:nSnpPerChr,10),
                  geno)

# Create "pheno" data.frame for rrBLUP
pheno = data.frame(gid=names(geno)[-(1:3)],
                   trait1=pop@pheno[,1])

# Get QTL names
qtl = getQtlMap()$id

## GWAS using kinship (K) adjustment
# Might take a while to run
modelK = GWAS(pheno, geno, plot=FALSE)

# Plot
manhattan(modelK, chr="chr", bp="pos", p="trait1", snp="snp",
          highlight=qtl, logp=FALSE, main="Adjusted for K")

@zjwinn above is an old GWAS script I have sitting around. I cleaned it up a little bit to be clearer. It doesn't address your original question, but I thought I'd share it in case it was useful.