Проход по столбцам в объектах S4 в R

Question

Я пытаюсь установить связь с помощью пакета snpStats .

У меня есть snp матрица под названием 'plink', которая содержит мои данные генотипа (в виде списка $ genotypes, $ map, $ fam), и plink $ genotype имеет: имена 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

У меня также есть фенотип данных, который содержит результаты (исход1, результат2, ...), которые я хотел бы связать сгенотип, такой:

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)

MyВопрос в том, как мне просмотреть результаты?Этот тип данных, кажется, отличается от всех других, и у меня возникают проблемы с манипулированием им, поэтому я также был бы признателен, если у вас есть предложения по некоторым учебникам, которые могут помочь мне понять, как сделать это, и другие манипуляции, такие как поднабор 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 - Дэвид Клэйтон.Хотя веб-сайт, указанный в описании пакета, неверен, он все еще находится в этом домене, и можно выполнить поиск документации с помощью функции расширенного поиска Google с этой спецификацией:

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

Вероятная причинаВаша проблема с доступом заключается в том, что это пакет S4 и методы доступа различны.Вместо методов печати у объектов S4 обычно есть show-методы.На упаковке есть виньетка: https://www -gene.cimr.cam.ac.uk / staff / clayton / курсы / 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

В документации сказано, что 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) или векторный фенотип (если он задан как простой вектор - т.е. в родительском кадре или в вашем «текущем» кадре, поскольку подпрограмма оценивается в «дочернем» кадре!)

Answer 3

Правильный способ сделать то, что требуется для первоначального плаката:

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) 
}