循环遍历R中的S4对象中的列

Question

我试图使用snpStats包执行关联。我有一个名为'plink'的snp矩阵，其中包含我的基因型数据（如 $基因型列表，$ map，$ fam），并且plink $ genotype具有：SNP名称作为列名称（2个SNP）和主题标识符作为行名称：

plink$genotype 
SnpMatrix with 6 rows and 2 columns 
Row names: 1 ... 6 
Col names: 203 204 

的砰砰数据集可以再现复制下面的PED和映射文件并将其保存为“plink.ped”分别plink.map”：

plink.ped: 

1 1 0 0 1 -9 A A G G 
2 2 0 0 2 -9 G A G G 
3 3 0 0 1 -9 A A G G 
4 4 0 0 1 -9 A A G G 
5 5 0 0 1 -9 A A G G 
6 6 0 0 2 -9 G A G G 

plink.map: 

1 203 0 792429 
2 204 0 819185 

然后以这种方式使用plink：

./plink --file plink --make-bed 

@----------------------------------------------------------@ 
|  PLINK!  |  v1.07  | 10/Aug/2009  | 
|----------------------------------------------------------| 
| (C) 2009 Shaun Purcell, GNU General Public License, v2 | 
|----------------------------------------------------------| 
| For documentation, citation & bug-report instructions: | 
|  http://pngu.mgh.harvard.edu/purcell/plink/  | 
@----------------------------------------------------------@ 

Web-based version check (--noweb to skip) 
Recent cached web-check found...Problem connecting to web 

Writing this text to log file [ plink.log ] 
Analysis started: Tue Nov 29 18:08:18 2011 

Options in effect: 
--file /ugi/home/claudiagiambartolomei/Desktop/plink 
--make-bed 

2 (of 2) markers to be included from [ /ugi/home/claudiagiambartolomei/Desktop /plink.map ] 
6 individuals read from [ /ugi/home/claudiagiambartolomei/Desktop/plink.ped ] 
0 individuals with nonmissing phenotypes 
Assuming a disease phenotype (1=unaff, 2=aff, 0=miss) 
Missing phenotype value is also -9 
0 cases, 0 controls and 6 missing 
4 males, 2 females, and 0 of unspecified sex 
Before frequency and genotyping pruning, there are 2 SNPs 
6 founders and 0 non-founders found 
Total genotyping rate in remaining individuals is 1 
0 SNPs failed missingness test (GENO > 1) 
0 SNPs failed frequency test (MAF < 0) 
After frequency and genotyping pruning, there are 2 SNPs 
After filtering, 0 cases, 0 controls and 6 missing 
After filtering, 4 males, 2 females, and 0 of unspecified sex 
Writing pedigree information to [ plink.fam ] 
Writing map (extended format) information to [ plink.bim ] 
Writing genotype bitfile to [ plink.bed ] 
Using (default) SNP-major mode 

Analysis finished: Tue Nov 29 18:08:18 2011 

我也有一个包含结果的表型数据帧（outcome1，outcome2，...）我想与基因型，这是该关联：

ID<- 1:6 
sex<- rep(1,6) 
age<- c(59,60,54,48,46,50) 
bmi<- c(26,28,22,20,23, NA) 
ldl<- c(5, 3, 5, 4, 2, NA) 
pheno<- data.frame(ID,sex,age,bmi,ldl) 

协会工程当我做的唯一条件是：（用公式“snp.rhs.test”）：

bmi<-snp.rhs.tests(bmi~sex+age,family="gaussian", data=pheno, snp.data=plink$genotype) 

我的问题是，我怎么通过成果循环？这种类型的数据 似乎不同于其他所有人，我在操作它时遇到问题 ，所以如果您有一些教程的建议 可以帮助我了解如何执行此操作以及其他操作（如子集），我也将不胜感激。例如，snp.matrix数据。

这是我曾尝试循环：

rhs <- function(x) { 
x<- snp.rhs.tests(x, family="gaussian", data=pheno, 
snp.data=plink$genotype) 
} 
res_ <- apply(pheno,2,rhs) 

Error in x$terms : $ operator is invalid for atomic vectors

然后我尝试这样的：

for (cov in names(pheno)) { 
association<-snp.rhs.tests(cov, family="gaussian",data=pheno, snp.data=plink$genotype) 
} 

Error in eval(expr, envir, enclos) : object 'bmi' not found

感谢您像往常一样对你的帮助！ -f

Answer 1

snpStats的作者是David Clayton。虽然在包装说明中列出的网站是错误的，他仍然在该领域，并有可能做一个搜索的文档与此规范的谷歌高级搜索功能：

snpStats site:https://www-gene.cimr.cam.ac.uk/staff/clayton/ 

可能的原因为您的难度与访问是这是一个S4包，访问的方法是不同的。而不是打印方法S4对象通常具有show-methods。这里有一个小插件：https://www-gene.cimr.cam.ac.uk/staff/clayton/courses/florence11/practicals/practical6.pdf，他的整个短期课程目录是开放的访问：https://www-gene.cimr.cam.ac.uk/staff/clayton/courses/florence11/

很明显，从snp.rhs.tests返回的对象可以用“[”使用序列号或名称，如第7页所示。你可以得到名称：

# Using the example on the help(snp.rhs.tests) page: 

> names(slt3) 
[1] "173760" "173761" "173762" "173767" "173769" "173770" "173772" "173774" 
[9] "173775" "173776" 

您可致电列中的东西可能是“插槽”

> getSlots(class(slt3)) 
    snp.names var.names  chisq   df   N 
     "ANY" "character" "numeric" "integer" "integer" 
> str(getSlots(class(slt3))) 
Named chr [1:5] "ANY" "character" "numeric" "integer" "integer" 
- attr(*, "names")= chr [1:5] "snp.names" "var.names" "chisq" "df" ... 
> names(getSlots(class(slt3))) 
[1] "snp.names" "var.names" "chisq"  "df"  "N"   

但对于遍历这些插槽名称没有[i,j]方法。您应该转到帮助页?"GlmTests-class"，其中列出了为该S4类定义的方法。

Answer 2

做的正确的方式所需的初始海报是什么：

for (i in ncol(pheno)) { 
    association <- snp.rhs.tests(pheno[,i], family="gaussian", snp.data=plink$genotype) 
} 

的snp.rhs.tests()的文件说，如果数据丢失，表型是从父帧拍摄 - 或许它是在措辞相反意义：如果指定了数据，表型将在指定的data.frame中进行评估。

这是一个清晰的版本：

for (i in ncol(pheno)) { 
    cc <- pheno[,i] 
    association <- snp.rhs.tests(cc, family="gaussian", snp.data=plink$genotype) 
} 

Answer 3

文档说data=parent.frame()是snp.rhs.tests()默认。

apply()代码有一个明显的错误 - 请不要做x <- some.fun(x)，因为它确实很糟糕。试试这个 - 删除data=，并使用不同的变量名称。

rhs <- function(x) { 
y<- snp.rhs.tests(x, family="gaussian", snp.data=plink$genotype) 
} 
res_ <- apply(pheno,2,rhs) 

此外，最初的海报的问题是误导。

plink $ genotype是一个S4对象，pheno是一个data.frame（一个S3对象）。您真的只想选择S3 data.frame中的列，但您将如何从snp.rhs.tests（）查找列（如果给出data.frame）或向量表型（if它是作为一个普通的矢量给出的 - 即在父框架或你的“当前”框架中，因为子程序是在“子框架”中评估的！）