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

Question

我正在尝试使用snpStats包执行关联。

我有一个名为'plink‘的snp矩阵，其中包含我的基因型数据(作为$genotypes、$map、$fam的列表)，并且plink$ subject有: SNP名称作为列名(2个SNP)，主题标识符作为行名：

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

可以复制plink数据集复制以下ped和map文件，并将它们分别保存为'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。尽管包描述中列出的网站是错误的，但他仍然在该域中，并且可以使用以下规范使用Google的高级搜索功能搜索文档：

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

访问困难的可能原因是这是一个S4包，访问方法是不同的。S4对象通常有显示方法，而不是打印方法。在这个包上有一个小插曲：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 )或向量表现型(如果它是以普通向量形式给出的--即在父框架或“当前”框架中，因为子例程是在“子”框架中求值的！)