phecodeX_PheWAS_example_R_script.txt

library(PheWAS)
set.seed(1)

## make sure you add the right path!
phecodeX_labels = read.csv('phecodeX_R_labels.csv')
phecodeX_rollup_map = read.csv('phecodeX_R_rollup_map.csv')
phecodeX_map = read.csv('phecodeX_R_map.csv') ## if you are using ICD-10 (not CM), load phecodeX_R_map_ICD_10_WHO.csv instead
phecodeX_sex = read.csv('phecodeX_R_sex.csv')

#Generate some example data
ex=generateExample()

#Extract the three parts from the returned list
id.vocab.code.count=ex$id.vocab.code.count
genotypes=ex$genotypes
id.sex=ex$id.sex

#Create the phecode table- translates the codes, adds
#exclusions, and reshapes to a wide format.
#Sum up the counts in the data where applicable.
#Note the vocabulary.map and rollup.map are set to
#phecodeX files
phenotypesX = createPhenotypes(id.vocab.code.count,
                               aggregate.fun=sum, id.sex=id.sex,
                               vocabulary.map=phecodeX_map,
                               rollup.map=phecodeX_rollup_map, 
                               sex.restriction=phecodeX_sex)

#Combine the data
data = inner_join(inner_join(phenotypesX,genotypes),id.sex)

#Run the PheWAS
results = phewas_ext(data,phenotypes=names(phenotypesX)[-1], genotypes=c("rsEXAMPLE"), covariates=c("sex"))

## Add phecode_labels
results_annotated = merge(phecodeX_labels, results, by='phenotype')

## Make a plot
phenotypeManhattan(results_annotated, 
                   title="My Example PheWAS Manhattan Plot",
                   annotate.phenotype.description=T,
                   sort.by.category.value=F, x.group.labels=T,
                   use.color=T)