[SOLVED] Problem using 'by' argument of 'descriptives.trait'

Questions about GenABEL (aka *ABEL) suite of packages
Forum rules
Please remember not to post any sensitive data on this public forum.
The first few posts of newly registered users will be moderated in order to filter out any spammers.

When get a solution to the problem you posted, please change the topic name (e.g. from "how to ..." to "[SOLVED] how to ..."). This will make it easier for the community to follow the posts yet to be attended.
richard.ahn
Bug reporter
Bug reporter
Posts: 8
Joined: Mon Jan 24, 2011 2:12 am

[SOLVED] Problem using 'by' argument of 'descriptives.trait'

Postby richard.ahn » Mon Jan 24, 2011 2:21 am

I'm using GenABEL 1.6-4 with R 2.12 on Ubuntu 10.04 to perform a GWAS with ~545K markers and 2258 samples, a binary disease variable, and age and sex as my other covariates. I'm running into a couple of issues which I can't seem to resolve and I thought someone on the board may be able to help me out.

First, in my phdata, I have a disease variable, bt, which is 0 and 1 only, no missing values (i.e., '.', ' ', or 'na') and no other values. Yet, when I run descriptives.trait with the 'by = bt' argument, I get the following:

Code: Select all

> descriptives.trait(data=cidr.pheno, by = bt)
Error in t.test.formula(ctrao ~ svar) :
  grouping factor must have exactly 2 levels
#I ran table to check the levels of the bt variable:
> table(cidr.pheno$bt)

   0    1
 530 1728


I also get the following message when I run qtscore with bt:

Code: Select all

Warning message:
In test.type(y, trait.type) : binomial trait is analysed as gaussian


The sum is 2258 which is the number of samples I have in my initial dataset. I verified that this was the case (n_bt = 2258) with other tools. Since the original phenotype file was a csv file from a windows machine, I also checked to make sure there were no dos characters in the file. Is there something that I need to do to transform the variable? I've gone through the GenABEL manual and tutorial and have also scoured the net for information and I'm not finding much success in resolving this issue.

Thank you!
Last edited by richard.ahn on Thu Jan 27, 2011 4:02 pm, edited 1 time in total.
Richard Ahn, M.A.
Graduate Student Researcher
Ph.D. Student
Department of Epidemiology
University of California, Irvine

yurii
GenABEL developer
GenABEL developer
Posts: 263
Joined: Fri Jan 21, 2011 5:20 pm

Re: GenABEL: Binary Phenotypic Variable Problem

Postby yurii » Mon Jan 24, 2011 3:00 pm

It is slightly difficult to tell what might went wrong here -- in general it would be great if you could send us a small example by running which we can reproduce an error.

Anyways, let us try to figure out what is going on and see if we can reproduce the error (and hopefully have an easy fix). Can you please run and send output for the following commands? I assume that class(cidr.pheno) is "gwaa.data"

class(cidr.pheno)
nids(cidr.pheno)
class(bt)
length(bt)
table(bt,exclude=NULL)
class(phdata(cidr.pheno)$bt)
table(phdata(cidr.pheno)$bt,exclude=NULL)

Yurii
Note that (Gen)ABELs are dynamically developing; while this post is intended to provide full information at the time of posting, please read on further posts, if any, as the topic may be updated with novel solutions at a later stage.

best regards,
Yurii

richard.ahn
Bug reporter
Bug reporter
Posts: 8
Joined: Mon Jan 24, 2011 2:12 am

Re: GenABEL: Binary Phenotypic Variable Problem

Postby richard.ahn » Mon Jan 24, 2011 4:17 pm

Code: Select all

> cidr.pheno = phdata(cidr)
> attach(cidr.pheno)
> class(cidr)
[1] "gwaa.data"
attr(,"package")
[1] "GenABEL"
> nids(cidr)
[1] 2258
> class(bt)
[1] "integer"
> length(bt)
[1] 2258
> table(bt,exclude=NULL)
bt
   0    1 <NA>
 530 1728    0
> class(cidr.pheno$bt)
[1] "integer"
> table(cidr.pheno$bt,exclude=NULL)

   0    1 <NA>
 530 1728    0


Do you think this problem arose because bt is an integer?
Richard Ahn, M.A.
Graduate Student Researcher
Ph.D. Student
Department of Epidemiology
University of California, Irvine

yurii
GenABEL developer
GenABEL developer
Posts: 263
Joined: Fri Jan 21, 2011 5:20 pm

Re: GenABEL: Binary Phenotypic Variable Problem

Postby yurii » Wed Jan 26, 2011 11:18 am

This is awkward. This is definitely not because it is integer, and I can not reproduce the error. See below the script which works just fine:

Code: Select all

library(GenABEL)
data(srdta)
xx <- 1*(runif(nids(srdta))>.5)
table(xx,exclude=NULL)
class(xx)
descriptives.trait(srdta,by=xx)
xx <- as.integer(xx)
class(xx)
descriptives.trait(srdta,by=xx)

What I suggest -- can you work out a minimal example (e.g. truncate your data, remove any sensitive information) and post the script which fails + the data it fails? The point here is that for the moment we can not reproduce the error and thus have no clue how to fix! If you do that, we then will put that in our bug tracker (https://r-forge.r-project.org/tracker/?atid=2058&group_id=505&func=browse)

In general -- this procedure is buggy in that it does not work always, sorry about that. We consider this procedure of relatively low importance, so it is 'on hold' still with few bugs reported already. May be your bug report will be the one after which someone will finally fix the procedure!

best,
Yurii
Note that (Gen)ABELs are dynamically developing; while this post is intended to provide full information at the time of posting, please read on further posts, if any, as the topic may be updated with novel solutions at a later stage.

best regards,
Yurii

richard.ahn
Bug reporter
Bug reporter
Posts: 8
Joined: Mon Jan 24, 2011 2:12 am

Re: GenABEL: Binary Phenotypic Variable Problem

Postby richard.ahn » Thu Jan 27, 2011 3:12 am

I will go ahead and truncate the data and post (it's all de-identified).

For the data, I will go ahead and post an illumina preferred marker file with 10 samples and 10 markers as well as the corresponding phenotype file. Will this do?

Thank you.
Richard Ahn, M.A.
Graduate Student Researcher
Ph.D. Student
Department of Epidemiology
University of California, Irvine

yurii
GenABEL developer
GenABEL developer
Posts: 263
Joined: Fri Jan 21, 2011 5:20 pm

Re: GenABEL: Binary Phenotypic Variable Problem

Postby yurii » Thu Jan 27, 2011 7:01 am

The failure occurs in

Code: Select all

> descriptives.trait(data=cidr.pheno, by = bt)
Error in t.test.formula(ctrao ~ svar) :
  grouping factor must have exactly 2 levels

and the cidr.pheno = phdata(cidr).

This means that there is no need to send genotypic data, only the minimal anonymized pheno-frame (like cidr.pheno) should do. Send it as an RData file (or a text file) + code which loads the data and brings us to the point when the 'descriptives.trait' command fails to do what it should do.

Yurii
Note that (Gen)ABELs are dynamically developing; while this post is intended to provide full information at the time of posting, please read on further posts, if any, as the topic may be updated with novel solutions at a later stage.

best regards,
Yurii

richard.ahn
Bug reporter
Bug reporter
Posts: 8
Joined: Mon Jan 24, 2011 2:12 am

Re: Problem using 'by' argument of 'descriptives.trait'

Postby richard.ahn » Fri Jan 28, 2011 2:16 am

Okay, here's my data:
The marker file: 'test_markers' (10 markers)
The phenotype file: 'test_phenos' (11 samples)
(They are located in test.zip)

Code: Select all

convert.snp.illumina(infile="test_markers", outfile="test.raw")
test = load.gwaa.data(pheno="test_phenos", geno="test.raw", force=T)
attach(phdata(test))
descriptives.trait(data=test, by = bt)


Please let me know if you have any problems with the data--I tested it before attaching.

*Yurii, I'm sorry I'm somehow missed your last post about not submitting the genotypes. The files are very small so I guess there's no harm done.

Richard
Attachments
test.zip
(1.05 KiB) Downloaded 417 times
Richard Ahn, M.A.
Graduate Student Researcher
Ph.D. Student
Department of Epidemiology
University of California, Irvine

yurii
GenABEL developer
GenABEL developer
Posts: 263
Joined: Fri Jan 21, 2011 5:20 pm

Re: Problem using 'by' argument of 'descriptives.trait'

Postby yurii » Fri Jan 28, 2011 7:01 am

Many thanks, it perfectly reproduces! I have added it to our bug tracker (with subj "[#1259] descriptives.trait with by.var argument does not work (contributed by Richard Ahn)"):

https://r-forge.r-project.org/tracker/i ... &atid=2058

I hope some volunteer contributor could fix that within soon -- we already have three bugs for this procedure, it is a leader!

Yurii
Note that (Gen)ABELs are dynamically developing; while this post is intended to provide full information at the time of posting, please read on further posts, if any, as the topic may be updated with novel solutions at a later stage.

best regards,
Yurii

richard.ahn
Bug reporter
Bug reporter
Posts: 8
Joined: Mon Jan 24, 2011 2:12 am

Re: Problem using 'by' argument of 'descriptives.trait'

Postby richard.ahn » Fri Jan 28, 2011 8:32 pm

Great! Hopefully, this bug will get resolved :)
Richard Ahn, M.A.
Graduate Student Researcher
Ph.D. Student
Department of Epidemiology
University of California, Irvine

Nicola Pirastu
GenABEL senior expert
GenABEL senior expert
Posts: 151
Joined: Wed Feb 09, 2011 3:24 pm

Re: Problem using 'by' argument of 'descriptives.trait'

Postby Nicola Pirastu » Mon Mar 07, 2011 11:07 pm

Hi, I figured out what the problem is:

basically all your controls are missing age so the t.test crushes, if this is all the data you're using that's the problem. When it says the 2 levels are required it is actually saying that you got only one level.
You can stil get results for the other colmuns with:

Code: Select all

descriptives.trait(data=test@phdata[,-3], by = (test@phdata$bt))


Do remember that when you use by.var you should put in the vector and not just the name of the variable or it won't work.

Hope I was of help.

Nicola


Return to “GenABEL”

Who is online

Users browsing this forum: No registered users and 2 guests