I am working on a genome-wide association study. I wish to examine the time to onset of a disease in relation to the alleles present at each SNP (survival analysis). I have done this for a few snps using STATA, but this won't work for a whole genome dataset so I am trying to do it in R.

My files are in PLINK format, and I have included the time to event data in the fam section of the read in plink data using the merge command.

My script is below:

```
#Read in plink files
GWAS<-read.plink("R.bed","R.bim","R.fam")
pheno <- read.delim("C:/Users/USER/Desktop/pheno.txt")
view(pheno)

library("GenABEL", lib.loc="~/R/win-library/3.1")

#find columns 
head(GWAS)

#merge plink file with followup and presence of disease
GWAS$fam <-merge(GWAS$fam,pheno, by.x="member", by.y= "Chip")

#this is now included in the plink files
#running in cox model using GenABEL package and mlreg command
coxm <- mlreg(GWAS((fam$followup),(map$allele.1)),"fam$disease")
```

However, it's not working, I'm pretty new at this and I'm not sure how to specify the genotype at each snp for each individual. Has anyone ran a genome-wide cox model in R? If so help would be greatly appreciated!
I figured it out, thank's everyone for the help. Thought I'd post up the script I used in case anyone else has similar issues.

What my problem was, was that I had the data in the wrong format, so I converted my bed, bim and fam files to ped and map using `--recode` in plink. The files could then be read into R and converted again to the gwaa.data class. Once it was read in the rest was straight forward.

```
#Load package
source("http://bioconductor.org/biocLite.R")
biocLite('GenABEL')

#open package
library('GenABEL')

#set working directory
setwd("C:/where_my_files_are")

#Convert to gwaa.data class
convert.snp.ped(ped='Rel2-0.ped',mapfile='Rel2-0.map', out='rel2-0.out')
survivaldata <- load.gwaa.data(phenofile = "pheno.txt", genofile = "rel2-0.out")

#run quality control
qc1 <- check.marker(survivaldata)
qcGWAS <- survivaldata[qc1$idok, qc1$snpok]

#Run cox proportional hazards model
coxresults <- mlreg(GASurv(survivaldata@phdata$followuptime, survivaldata@phdata$event)~survivaldata@phdata$age, data=qcGWAS)

#View top 20 results
descriptives.scan(results, top=20)
```