Sex as an Interactive Covariate

[Adrianna Jurek]

Hello!

I have been using A Guide to QTL mapping to help me run sex as an additive and interactive covariate, and I had a couple questions relating to the code and the output. 

In the book, sex and cross are used as covariates, but for my data, I only have sex as a covariate. I modified the code as seen below, but I wasn't sure if I handled the interactive covariate correctly. When I run the interactive covariate code to determine LOD significance, the LOD score is 151 for the 5% threshold. I'm not sure if that is the correct score and I'm overthinking it, or if I handled the code incorrectly leading to an inflated LOD threshold. 

Here is the code I used:
#addcovar
#use mydataf1 for step.genoprob
sex<-mydataf1$pheno$sex
out.0<-scanone(mydataf1, method="hk")
out.a<-scanone(mydataf1, addcovar=sex, method="hk")
plot(out.a,out.0, ylab="LOD score", col=c("red","black"),alternate.chrid=TRUE,bandcol="gray87")
operm.0<-scanone(mydataf1, n.perm=1000, perm.Xsp=TRUE)
operm.a<-scanone(mydataf1, addcovar=sex, n.perm=1000, perm.Xsp=TRUE)
out.both<-c(out.0,out.a, labels=c("nocovar","covar"))
operm.both<-cbind(operm.0, operm.a, labels=c("nocovar","covar"))

summary(out.both, perms=operm.both, format="allpeaks", alpha=0.5, pvalues=TRUE)

out.i<-scanone(mydataf1, addcovar=sex, intcovar=sex, method="hk")

plot(out.i, out.i-out.a, ylab="LOD Score", col=c("blue", "red"), alternate.chrid=TRUE)

set.seed(54955149)
operm.a<-scanone(mydataf1, addcovar=sex, n.perm=1000, perm.Xsp=TRUE, verbose=FALSE)
set.seed(54955149)
operm.i<-scanone(mydataf1, addcovar=sex, intcovar=sex, n.perm=1000, perm.Xsp=TRUE, verbose=FALSE) 

out.ia<-c(out.i, out.i-out.a, labels=c("f", "i"))
operm.ia<-cbind(operm.i, operm.i-operm.a, labels=c("f", "i"))

summary(out.ia, perms=operm.ia, alpha=0.05, pvalues=TRUE, lodcolumn=2)

Output for last line of code: 
Autosome LOD thresholds (1000 permutations)
   lod.f lod.i
5%   151   148

X chromosome LOD thresholds (14561 permutations)
   lod.f lod.i
5%  3.84  2.77

There is a clear sex difference in the data, but I want to make sure I am running and interpreting the results from the addcovar and intcovar runs correctly.

Thank you for your time!
Adrianna 

[Karl Broman]

The code looks fine, but yeah, the LOD thresholds of 151 and 148 are crazy high.

I would look at plots of the permutation results.

I'm also surprised that the X chromosome threshold for lod.i is not 0. 

Typically the X chromosome includes a sex as both an additive and interactive covariate, and so you get warnings like "Dropped 1 additive covariate on X chromosome" and "Dropped 1 interactive covariate on X chromosome." and then the results are the same on the X chromosome, whether you include sex as a covariate or not. And that should hold for all permutations.

karl

[Adrianna Jurek]

Hi Dr. Broman,

I re-ran the code for the thresholds and I did get the lod.i to be 0 this time, not sure where I went wrong with the first one but I have been editing the code to come up with the one I sent you so I may have run the wrong code line. Here are the updated values:  Autosome LOD thresholds (1000 permutations)

   lod.f lod.i
5%   151   148

X chromosome LOD thresholds (14561 permutations)
   lod.f lod.i

5%  3.83     0

As for looking at the permutation plots, the addcovar graph looks quite different from the intcovar graph with a bar at 150 and the rest in the normal range. These (attached) were graphed using plot(operm.i) and plot(operm.a). I'm not sure what to make of the difference here. Any insight would be appreciated. 

Thanks for your help, 

Adrianna

[Karl Broman]

More than 5% of the time, your permutations are giving you an absurdly large LOD score, when sex is included as an interactive covariate.

So the next task is to figure out why that's the case. You might try doing the permutations "by hand" one at a time, and look for one that's giving a huge LOD score, and then studying that case.

karl

[Adrianna Jurek]

Hi Dr. Broman, 

I ran the permutations "by hand" to generate a large lod score, and I got the results below. I did notice that this run in particular had more warnings than the other runs (warnings 1-6 do not show up in any other permutation that I ran). I'm not sure exactly how to analyze past this, though.  Do you have any recommendations?

> operm.i<-scanone(mydataf1, intcovar=sex,addcovar=sex, n.perm=1, perm.Xsp=TRUE, verbose=FALSE)
Warning messages:
1: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
2: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
3: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
4: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
5: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
6: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  X'X matrix is singular.
7: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  Dropped 1 additive covariates on X chromosome.
8: In scanone(cross, , pheno.col, model, method, addcovarp, intcovarp,  :
  Dropped 1 interactive covariates on X chromosome.
> operm.i
$A
     lod
1 151.39

$X
         lod
1  1.2136137
2  0.7784477
3  1.3997147
4  1.5140616
5  0.9615083
6  1.1959332
7  2.6498336
8  0.5812863
9  1.0904050
10 1.1189155
11 0.4472567
12 0.5589734
13 0.4677648
14 0.8206034
15 2.1426072

attr(,"xchr")
    1     2     3     4     5     6     7     8     9    10    11    12
FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
   13    14    15    16    17    18    19     X
FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
attr(,"L")
         A          X
1185.68900   81.43041
attr(,"method")
[1] "em"
attr(,"model")
[1] "normal"
attr(,"type")
[1] "f2"
attr(,"class")
[1] "scanoneperm" "list"      
> summary(operm.i)
Autosome LOD thresholds (1 permutations)
    lod
5%  151
10% 151

X chromosome LOD thresholds (15 permutations)
     lod
5%  2.63
10% 2.60

Thanks,
Adrianna

[Karl Broman]

What I meant, in doing it by hand, is something like the following:

for(i in 1:1000) {
    tempx <- mydataf1
    tempx$pheno <- mydataf1$pheno[sample(nrow(mydataf1$pheno)),,drop=FALSE]
    tempsex <- tempx$pheno$sex
    out_perm <- scanone(tempx, method="hk", addcovar=tempsex, intcovar=tempsex)
    m <- max(out_perm[,3])
    cat(i, m, "\n")
    if(m > 100) break
}

When it stops, tempx would be a cross with the phenotype rows in the order that produced the LOD>100 with sex as an interactive covariate.

karl