Meta-analytic Mendelian randomization (MR) of expression of cytosolic and mitochondrial RPs on breast cancer (BCAC 2017) in blood
MR-dysreg-GTEx-all-RP-BC-blood
2021-04-18
Abstract
Under conditions of nucleolar stress, various ribosomal proteins (RPs) interact with the MDM2 proto-oncogene (MDM2) to stabilize p53. Dysregulation of the RPs may disturb this pathway, but whether changes in the expression of the RPs influence risk for breast cancer in the general population is unknown. We tested whether expression changes in cytosolic or mitochondrial RPs influences risk for breast cancer using inverse-variance weighted (IVW) Mendelian randomization (MR)–an instrumental variables technique. Summary statistics for expression quantitative trait loci (eQTLs) for SNPs decreasing expression of RPs (P<5x10-06) served as instrumental variables for decreased RP expression. The summary statistics for the SNPs decreasing expression of the RPs were extracted from the Gene-Tissue Expression (GTEx) project’s (version 8 “eGene”) eQTL data in blood tissue. Summary statistics for these same SNPs were extracted from the Breast Cancer Association Consortium’s (BCAC) 2017 genome-wide association study (GWAS) of breast cancer (122977 cases and 105974 controls). After outliers were removed, summary statistics for 34 SNPs were used to ascertain whether changes in expression of RPs increases risk for breast cancer. This was done by calculating and meta-analyzing individual Wald ratios. Wald ratios were obtained by dividing the effect estimate for the impact of a given SNP on breast cancer by the effect estimate for the impact of the same SNP on RP expression. The meta-analyzed Wald ratios yield an overall causal (IVW) estimate that should be free from most sources of environmental confounding, so long as the assumptions for MR are not violated. Chief among the MR assumptions are that the genetic instruments (RP eQTLs) influence breast cancer only through the expression of the RPs and not through pleiotropic pathways. We performed sensitivity analyses, using the MR-Egger, weighted median, and weighted mode MR estimators to test for pleiotropy. The IVW estimate revealed an increased risk for breast cancer corresponding to an increase in expression of cytosolic and mitochondrial RPs, regardless of whether the effect alleles increased or decreased expression (odds ratio=1.011; 95% CI=1.000, 1.021; P=0.042). There was no evidence of pleiotropy. The results implies that environmental factors that increase the expression of the RPs in blood may increase risk for breast cancer and/or changes in RP expression in blood act as biomarkers of RP expression changes in breast. However, the confidence interval essentially includes 1, even though the p-value is <0.05. This means that the finding is very close to null.
Contents
0.1 Overview
- Method: meta-analytic Mendelian randomization (MR)
- Instrument data source: Gene-Tissue Expression (GTex) Project
- Expression quantitative trait loci (eQTL)
- Outcome data source: BCAC 2017 release
- Ribosomal protein (RP): small & large subunits only
- Significance threshold for analyses here (Q-value for ‘eGene’ status <0.05 & nominal P<5x10-6)
0.2 Actions
# Create the directory to match the repository in GitHub
mkdir /n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood
cd /n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood- Read in the full GTEx v8 “eGene” data
- Clean up column names
- Extract blood eQTLs
- Extract RP genes
- Select SNPs at P<5e-06
- Format for MR
- Merge blood RP eqtls with BCAC
- Run MR
- Look for outliers & remove them
- Test for pleiotropy & heterogeneity
- Generate forest plot & scatterplot
- Table the data sources
0.3 Code
# MR of any RP expression changes in blood on overall breast cancer
setwd('/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood')
# Hand-curated RP list (no pseudogenes): 150
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
rp_list_expanded
# Look and clean; highlighted out but very important
# master_v8_n=read.csv("/n/home04/cdadams/ld/ldsc/v8_master_tissues_n.csv")
# head(master_v8_n)
# master_v8_n$X=NULL
# colnames(master_v8_n)[2] <- "chromosome"
# colnames(master_v8_n)[31] <- "MAF"
# master_v8_n$eaf=ifelse(master_v8_n$ref_factor=='1', master_v8_n$MAF, 1-master_v8_n$MAF)
# gtex_blood=master_v8_n[which(master_v8_n$tissue=="Whole_Blood"),]
# write.csv(gtex_blood, '/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/gtex_blood.csv')
gtex_blood=read.csv('/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/gtex_blood.csv')
gtex_blood$X=NULL
str(gtex_blood)
str(rp_list_expanded)
rp_list_expanded$all_RP=as.factor(rp_list_expanded$all_RP)
gtex_blood$gene_name=as.factor(gtex_blood$gene_name )
gtex_blood_RP=gtex_blood[which(gtex_blood$gene_name %in% rp_list_expanded$all_RP),]
dim(gtex_blood_RP)
gtex_blood_RPClean <- gtex_blood_RP[!duplicated(gtex_blood_RP$gene_name),]
dim(gtex_blood_RPClean)
table(gtex_blood_RPClean$SNP)
gtex_blood_RPClean$gene_name=droplevels(gtex_blood_RPClean$gene_name)
freq1=table(gtex_blood_RPClean$gene_name)
freq1=as.data.frame(freq1)
freq1 #150!
dim(freq1)
sum(freq1$Freq) #150
# Get the RP eqtls
# Select the P<5x10-06 eqtls
gtex_blood_RPClean_p_06=gtex_blood_RPClean[which(gtex_blood_RPClean$pval<0.000005),]
freq2=table(gtex_blood_RPClean_p_06$gene_name)
freq2=as.data.frame(freq2)
freq2
sum(freq2$Freq) #69
# # Order
gtex_blood_RPClean_p_06=gtex_blood_RPClean_p_06[order(gtex_blood_RPClean_p_06$gene_name, decreasing=TRUE),]
write.csv(gtex_blood_RPClean_p_06, '/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/gtex_blood_RPClean_p_06.csv')
# MR format
# Use samplesize_col = 'position' to keep position. Change after.
exposure_dat <- read_exposure_data(
filename = "/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/gtex_blood_RPClean_p_06.csv",
sep = ',',
snp_col = 'SNP',
beta_col = 'beta',
se_col = 'se',
effect_allele_col = 'effect_allele',
phenotype_col = 'tissue',
other_allele_col = 'other_allele',
eaf_col = 'eaf',
samplesize_col = 'position',
gene_col = 'gene_name',
pval_col = 'pval',
)
freq3=table(exposure_dat$gene.exposure)
freq3=as.data.frame(freq3)
dim(freq3)
sum(freq3$Freq)
freq3
# Format the colname names
exposure_dat$chr_name=exposure_dat$chr.exposure
exposure_dat$chrom_start=exposure_dat$samplesize.exposure
exposure_dat$samplesize.exposure=NULL
exposure_dat$gene_tissue=exposure_dat$exposure
exposure_dat$exposure="RP expression in blood"
exposure_dat$id.exposure="RP expression in blood"
# # Select the decreased expression eQTLs
# dat_decrease=exposure_dat[which(exposure_dat$beta.exposure<0),]
#
# # Subset by just the cytosolic RPs
# decrease_RPs <- dat_decrease[startsWith(as.character(dat_decrease$gene.exposure), 'R'),]
# Extract the increased expression loci from BCAC
outcome_dat <- extract_outcome_data(exposure_dat$SNP, c('1126'),
proxies = 1, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3)
dat <- harmonise_data(exposure_dat, outcome_dat, action = 2)
freq4=table(dat$gene.exposure)
freq4=as.data.frame(freq4)
freq4
sum(freq4$Freq)
dim(freq4)
# Clump for LD
dat_clump <- clump_data(dat)
freq5=table(dat_clump$gene.exposure)
freq5=as.data.frame(freq5)
freq5
sum(freq5$Freq)
dim(freq5)
# Run MR
mr_results<- mr(dat_clump)
mr_results
dat=dat_clump
# Detect and remove outliers (outliers can indicate pleiotropy)
raddat <- format_radial(BXG=dat$beta.exposure, BYG=dat$beta.outcome,
seBXG=dat$se.exposure, seBYG=dat$se.outcome, RSID=dat$SNP)
ivwrad <- ivw_radial(raddat, alpha=0.05, weights=1)
ivwrad2 <- ivw_radial(raddat, alpha=0.05, weights=1)
dim(ivwrad$outliers)[1]
outliers=ivwrad$outliers[1]
outliers
myvars=outliers$SNP
dat2 <- dat[! dat$SNP %in% myvars, ]
dim(dat2)
dim(dat)
dat2$outcome="breast cancer (BCAC 2017)"
dat2$id.outcome="breast cancer (BCAC 2017)"
# Run MR without outliers
mr_results <- mr(dat2)
mr_results
singlesnp_results=mr_singlesnp(dat2, parameters = default_parameters(),
single_method = "mr_wald_ratio",
all_method = c("mr_ivw",
"mr_egger_regression",
"mr_weighted_mode",
"mr_weighted_median"))
# Remove the snps that weren't used in the MR to have a record of the data source
singlesnp_results.2=singlesnp_results[1:34,]
dim(singlesnp_results.2)
dat2.1=dat2[which(dat2$SNP %in% singlesnp_results.2$SNP),]
dim(dat2.1)
dim(dat)
dat=dat2.1
dim(dat)
# r2 and Fstat
n=480
k=34
summary(dat$eaf.exposure)
test=scale(dat$beta.exposure)
str(test)
test <- scale(dat$beta.exposure)
test3=as.data.frame(test)
dim(test3)
dat$scaled.dat=test3$V1
summary(dat$scaled.dat)
p=dat$eaf.exposure
r2.1<-(2*(dat$scaled.dat^2)*p*(1-p))/1 #dat$beta.exposure
r2=sum(r2.1, na.rm = TRUE)
r2
#r2=mean(r2.1, na.rm=TRUE)
#r2
Fstat <- r2*(n-1-k)/((1-r2)*k)
Fstat
# Combine results
sin=singlesnp_results
res=mr_results
het<-mr_heterogeneity(dat)
plt<-mr_pleiotropy_test(dat)
# Bug in function for combining; save separately
# all_res<-combine_all_mrresults(res,het,plt,sin,ao_slc=F,Exp=T,split.exposure=T,split.outcome=T)
# all_res$lower_intercept=all_res$intercept-1.96*all_res$intercept_se
# all_res$upper_intercept=all_res$intercept+1.96*all_res$intercept_se
#
# head(all_res[,c("Method","outcome","exposure","nsnp","b","se","or", "or_lci95",
# "or_uci95","pval","intercept","intercept_se","intercept_pval",
# "Q","Q_df","Q_pval")])
#write.csv(all_res,"combined_results_all_RP_MRP_blood_v8.csv")
write.csv(mr_results,'/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/res_all_RP_blood_v8.csv')
write.csv(singlesnp_results,'/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/singlesnp_all_RP_blood_v8.csv')
write.csv(dat, '/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/all_RP_blood_v8_SIMEX.csv')
write.csv(plt,'/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/plt_all_RP_blood_v8.csv')
write.csv(het,'/n/home04/cdadams/Aim2/MR-dysreg-GTEx-all-RP-BC-blood/het_all_RP_blood_v8.csv')0.4 Figures
The top plot displays the individual SNPs included in the analysis, their effect estimates and 95% confidence intervals, along with the meta-analyzed “causal” estimates (the term “causal” is used in MR to mean MR is less prone to confounding and reverse causation than observational designs).
The inverse-variance weighted (IVW) estimate is the approach that is used to test the causal effect. The weighted median, weighted mode, and MR-Egger are sensitivity estimators. If the directions and magnitudes of the sensitivity estimators’ findings comport with those of the IVW, this is evidence against pleiotropy. (The sensitivity estimators needn’t all be statistically signficiant. For instance, the MR-Egger estimate isn’t, but this only means that the MR-Egger test doesn’t provide additional support for the IVW. The direction and magnitude of the MR-Egger estimate does, however, comport with the IVW’s.)
The bottom (scatter) plot is a sensitivity analysis for heterogeneity in the effect estimates. For the test of any RP expression on breast cancer, it shows a) the that SNPs changing the expression of the RPs increased risk for breast cancer across all estimators (the rising regression lines), b) that the sensitivity estimators comport (the regression lines don’t crisscross), and c) that there aren’t substantial outliers (there are no substantial outliers, since I detected and removed them before running the final MR; outliers can indicate pleiotropy).
The MR findings support an increased risk for breast cancer corresponding to increased expression of cytosolic and mitochondrial RPs in blood tissue. The IVW estimate’s confidence interval does not include 1.
FALSE Saving 8 x 5 in image
FALSE Saving 8 x 5 in image0.5 IVW & sensitivity estimator results
## method nsnp or se or_lci95 or_uci95 pval
## 1 MR Egger 34 1.001 0.010 0.982 1.021 0.886
## 2 Weighted median 34 1.002 0.008 0.987 1.018 0.753
## 3 Inverse variance weighted 34 1.011 0.005 1.000 1.021 0.042
## 4 Simple mode 34 0.998 0.015 0.969 1.027 0.875
## 5 Weighted mode 34 1.000 0.007 0.986 1.015 0.956