Wednesday, January 28, 2009

GenABEL

library(GenABEL)
data(srdta)

phdata<-srdta@phdata

library(genetics)
summary(datSNP$snp10001)
summary(datSNP$snp10002)
LD(datSNP$snp10001,datSNP$snp10002)
table(datSNP$snp10001,datSNP$snp10002)

### Structure of gwaa.data
# Phenotypic data
srdta@phdata
# number of SNPs
srdta@gtdata@nsnps
# IDs of SNPs
srdta@gtdata@snpnames
# Chromosome
srdta@gtdata@chromosome
# SNPs map positions
srdta@gtdata@map
# number of subjects
srdta@gtdata@nids
# IDs of subjects
srdta@gtdata@idnames
# gender
srdta@gtdata@male
# Genotypic data
srdta@gtdata@gtps

### exercise
# SNPs after 2,490,000 pb
srdta@gtdata@snpnames[srdta@gtdata@map>2490000]

### Sub-setting

# the first 5 subject and 3 SNPs
ssubs<-srdta[1:5,1:3]
srdta[c("p141","p147","p2000"),]@phdata
srdta[c("p141","p147","p2000"),c("rs10","rs29")]

# the first character is the consensus
as.character(ssubs@gtdata@coding)
# homozygotes for the first character is "0"
as.character(ssubs@gtdata@strand)

as.character(ssubs@gtdata)
as.numeric(ssubs@gtdata)

snps<-as.character(srdta@gtdata)

# convert data to the format for "genetics" library
as.genotype(ssubs@gtdata)

# convert data to the format for "haplo.stats" library
as.hsgeno(ssubs@gtdata)

# exercise 26
table(as.character(srdta[,"rs114"]))

### Exploring
summary(ssubs)