Plot of effect estimate comparisions as used in PsA immunochip paper

effect_estimate_comp.R

setwd("C:/Users/mdeasjdb/Dropbox/PsA_immunochip/effect_comparison")

# format known associtions table - this contains a link between the two results sets and r2
xref <- read.csv(file="psoriasis_known_loci.txt", header=TRUE)
colnames(xref)[6] <- 'xref_study_snp'
colnames(xref)[1] <- 'xref_tsoi_snp'
xref <- subset(xref, select=c('xref_study_snp', 'xref_tsoi_snp', 'ld_rs', 'gene'))
xref$ld_rs[is.na(xref$ld_rs)] <- 1
xref$ld <- as.numeric(as.character(xref$ld_rs))

# format Tsoi results to include confidence intervals.
tsoi <- read.table(file="tsoi_table.txt", header=TRUE)
tsoi <- tsoi[tsoi$snp != 'rs4406273',] # drop MHC snp
tsoi <- subset(tsoi, select=c('snp', 'p_value', 'or'))
colnames(tsoi) <- c('tsoi_snp', 'tsoi_p', 'tsoi_or')

# calculate confidence intervals
tsoi$coeff <- log(tsoi$tsoi_or)
tsoi$z <- -(qnorm(tsoi$tsoi_p/2))
tsoi$se <- tsoi$coeff/tsoi$z
tsoi$lci <- tsoi$coeff - (1.96*tsoi$se)
tsoi$uci <- tsoi$coeff + (1.96*tsoi$se)
tsoi$tsoi_l95 <- exp(tsoi$lci)
tsoi$tsoi_u95 <- exp(tsoi$uci)
tsoi <- subset(tsoi, select=c(colnames(tsoi)[grep('tsoi', colnames(tsoi))]))
               
# read study results
study <- read.table(file="known_loci_risk_allele_coded.assoc.logistic", header=TRUE)
study <- subset(study, select=c('SNP', 'OR', 'L95', 'U95', 'P'))
colnames(study) <- c('study_snp', 'study_or', 'study_l95', 'study_u95', 'study_p')

# merge data
data <- merge(study, xref, by.x='study_snp', by.y='xref_study_snp')
data <- merge(tsoi, data, by.x='tsoi_snp', by.y='xref_tsoi_snp')
data <- subset(data, ld >= 0.8)

# limit to assoicated in both?

# plot comparison
plot(data$tsoi_or, data$study_or, ylim=c(1,2), xlim=c(1,1.8), pch=19, col="darkgreen", main="Comparison of effect sizes between psoriasis and PsA", xlab="Psoriasis Odds Ratio (confidence intervals)", ylab="PsA Odds Ratio (confidence intervals)", log="xy", axes=FALSE)
axis(side=1)
axis(side=2)
grid()
abline(0,1, lty=2, col="darkred")
arrows(data$tsoi_l95, data$study_or, data$tsoi_u95, data$study_or, angle=90, code=3, length=0.05, col="grey", lwd=2)
arrows(data$tsoi_or, data$study_l95, data$tsoi_or, data$study_u95, angle=90, code=3, length=0.05, col="grey", lwd=2)
points(data$tsoi_or, data$study_or, pch=19, col="darkgreen")
text(1.35, 1.75, "IL23A")
text(1.3, 1.08, "LCE3B")