[GenABEL-dev] extracting results from OmicABEL runs - something is wrong

[Yurii Aulchenko]

Dear Sodbo, Diego,

It looks like either 'reshuffle' or CLACK-GWAS do something wrong - and I
actually think this is reshuffle (see my tests below). Basically, I run
analysis on 'big' genotype set (gtDataChr1) and a thinned one, where I take
every 100th marker only (gtDataChr1small). When I extract results with
reshuffle, I see a marker with chi2>25 in small, but not in the original
big data set!

You can find the test-files in my public dropbox folder (I am going to send
the link to you separately; if anyone else is interested in testing, let me
know); the commands I used for testing are provided below.

best, and let mw know,
YA

Running analysis:

CLAK-GWAS -cov covData -phi gKinId -pheno phDataChr1 -snp gtDataChr1small
-out BAllSmall -var eigen -nths 4
reshuffle BAllSmall --chi=25
[identifies one SNP-trait pair, rs3738547-GI_7657133.S, with chi2>25]
CLAK-GWAS -cov covData -phi gKinId -pheno phDataChr1 -snp gtDataChr1 -out
BAll -var eigen -nths 4
reshuffle BAll --chi=25
[no SNP-trait pairs identified at chi2>25]

To make sure that actually the data are right, I ran the following tests in
R making sure that the SNP is present in both big and thinned data set

> library(DatABEL)
DatABEL v.0.9-4 (March 12, 2013) loaded

> g <- databel("gtDataChr1")
> gs <- databel("gtDataChr1small")
> p <- databel("phDataChr1")
> snp <- as(g[,"rs3738547"],"vector")
> snps <- as(gs[,"rs3738547"],"vector")
> ph <- as(p[,"GI_7657133.S"],"vector")
> summary(lm(ph~snp))

Call:
lm(formula = ph ~ snp)

Residuals:
     Min       1Q   Median       3Q      Max
-0.33692 -0.04485  0.00333  0.05632  0.13978

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.009116   0.005337   1.708   0.0892 .
snp         -0.091529   0.015967  -5.732 3.53e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07357 on 204 degrees of freedom
Multiple R-squared:  0.1387, Adjusted R-squared:  0.1345
F-statistic: 32.86 on 1 and 204 DF,  p-value: 3.527e-08

> summary(lm(ph~snps))

Call:
lm(formula = ph ~ snps)

Residuals:
     Min       1Q   Median       3Q      Max
-0.33692 -0.04485  0.00333  0.05632  0.13978

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.009116   0.005337   1.708   0.0892 .
snps        -0.091529   0.015967  -5.732 3.53e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07357 on 204 degrees of freedom
Multiple R-squared:  0.1387, Adjusted R-squared:  0.1345
F-statistic: 32.86 on 1 and 204 DF,  p-value: 3.527e-08



And this is why I blame this on reshuffle:

yuryaulchenko$ perl extractCell.pl BAll 1029 115001
 0.00911574810743332 -0.0915286764502525 0.00573702435940504
0.0171629786491394 -2.74352405540412e-05
YAMacBookPro:analysisOfHapMapeQTL yuryaulchenko$ perl extractCell.pl
BAllSmall 1029 1151
 0.00911574810743332 -0.0915286764502525 0.00573702435940504
0.0171629786491394 -2.74352405540412e-05
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20130629/79a3f9ae/attachment.html>