Tuesday, February 3, 2009

haplo.stats

library(haplo.stats)

data(hla.demo); attach(hla.demo)

geno<-hla.demo[,c(17,18,21:24)]
label<-c("DQB","DRB","B")

geno.desc<-summaryGeno(geno, miss.val=c(0,NA))
print(geno.desc[c(1:10,80:85,135:140),])

table(geno.desc[,3])
miss.all<-which(geno.desc[,3]==3)
hla.demo.updated<-hla.demo[-miss.all,]

save.em<-haplo.em(geno=geno,locus.label=label,miss.val=c(0,NA))

summary(save.em, nlines=7)

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)