#cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR
#ssh cdadams@login.rc.fas.harvard.edu
#salloc -n 1 -p shared -t 0-09:00 --mem=9000–SMR results for complex diseases on a shiny app: https://shiny.cnsgenomics.com/SMRdb/
–Paper with helpful description of SMR: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4979185/
–Autism data come from: https://figshare.com/articles/dataset/asd2019/14671989?file=28169289
–This is the GWAS results file from the meta-analysis of ASD by the Lundbeck Foundation Initiative for Integrative Psychiatric Research (iPSYCH) and the Psychiatric Genomics Consortium (PGC) released in November 2017.
–Citation for studies using these data:
–Preprint by Grove et al: ‘Common risk variants identified in autism spectrum disorder’:
–https://www.biorxiv.org/content/early/2017/11/27/224774 (2017)
–Later published here:
–https://www.nature.com/articles/s41588-019-0344-8
–18,381 individuals with ASD and 27,969 controls
–Prepare the ASD vcf
–Get the ASD data
–Year 2017
–Category Disease
–Sub category Psychiatric / neurological
–Population European
–Sex Males and Females
–ncase 18,382
–ncontrol 27,969
–Sample size 46,351
–Number of SNPs 9,112,386
–Unit \N
–Author NA
–Consortium iPSYCH-PGC
–Ontology NA
–Build HG19/GRCh37
<!-- wget https://gwas.mrcieu.ac.uk/files/ieu-a-1185/ieu-a-1185.vcf.gz --># Read two times the vcf file, first for the columns names, second for the data
tmp_vcf<-readLines(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene/Brain_Amygdala/ieu-a-1185.vcf.gz"))
tmp_vcf_data<-read.table(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene/Brain_Amygdala/ieu-a-1185.vcf.gz"), stringsAsFactors = FALSE)
# filter for the columns names
tmp_vcf<-tmp_vcf[-(grep("#CHROM",tmp_vcf)+1):-(length(tmp_vcf))]
vcf_names<-unlist(strsplit(tmp_vcf[length(tmp_vcf)],"\t"))
names(tmp_vcf_data)<-vcf_names
head(tmp_vcf_data) #ES:SE:LP:AF:SS:ID
colnames(tmp_vcf_data)[1]="chromosome"
colnames(tmp_vcf_data)[2]="POS"
colnames(tmp_vcf_data)[3]="SNP"
colnames(tmp_vcf_data)[4]="A2"
colnames(tmp_vcf_data)[5]="A1"
colnames(tmp_vcf_data)[10]="mix"
tmp_vcf_data$b=sub("\\:.*", "", tmp_vcf_data$mix)
tmp_vcf_data$SE.1=str_before_nth(tmp_vcf_data$mix, ":", 2)
tmp_vcf_data$se=sub('.*:', '', tmp_vcf_data$SE.1)
tmp_vcf_data$p.1=str_before_nth(tmp_vcf_data$mix, ":", 3)
tmp_vcf_data$p=sub('.*:', '', tmp_vcf_data$p.1)
tmp_vcf_data$eaf.1=str_before_nth(tmp_vcf_data$mix, ":", 4)
tmp_vcf_data$eaf=sub('.*:', '', tmp_vcf_data$eaf.1)
# tmp_vcf_data$N.1=str_before_nth(tmp_vcf_data$mix, ":", 5)
# tmp_vcf_data$N=sub('.*:', '', tmp_vcf_data$N.1)
tmp_vcf_data$p=as.numeric(as.character(tmp_vcf_data$p))
tmp_vcf_data$b=as.numeric(as.character(tmp_vcf_data$b))
#tmp_vcf_data$N=as.numeric(as.character(tmp_vcf_data$N))
#tmp_vcf_data$eaf=as.numeric(as.character(tmp_vcf_data$eaf))
tmp_vcf_data$se=as.numeric(as.character(tmp_vcf_data$se))
summary(tmp_vcf_data$p)
tmp_vcf_data$pval=10^(-tmp_vcf_data$p)
tmp_vcf_data$pval2=10^(tmp_vcf_data$p)
summary(tmp_vcf_data$pval)
summary(tmp_vcf_data$p)
tmp_vcf_data$p=tmp_vcf_data$pval
tmp_vcf_data$n=46351
tmp_vcf_data$freq=NA
# Select asd variables for SMR
myvars=c("SNP", "A1", "A2", "freq", "b", "se", "p", "n" )
smr_asd=tmp_vcf_data[myvars]
smr_asd=smr_asd[order(smr_asd$p),]
write.table(smr_asd, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene/Brain_Amygdala/smr_asd.txt",
sep="\t", row.names = FALSE, col.names = TRUE)# Remove quotes to run in SMR
<!-- sed 's/\"//g' smr_asd.txt > smr_asd -->–README for cis-QTLs for the v8 GTEx release-
–This document describes the output format and column headers
–for cis-eQTL and cis-sQTL files from the v8 GTEx data release.
–eGenes
–File extension: .egenes.txt.gz, .sgenes.txt.gz –Column headers:
–gene_id: GENCODE/Ensembl gene ID
–gene_name: GENCODE gene name
–gene_chr: chromosome (gene)
–gene_start: gene start position (in base pairs; 1-based coordinates)
–gene_end: gene end position (in base pairs; 1-based coordinates)
–strand: genomic strand
–num_var: number of variants in cis-window
–beta_shape1: 1st shape parameter of the fitted Beta distribution: B(shape1, shape2)
–beta_shape2: 2nd shape parameter of the fitted Beta distribution: B(shape1, shape2)
–true_df: Effective degrees of freedom the Beta distribution approximation
–variant_id: variant ID in the format {chr}{pos_first_ref_base}{ref_seq}_{alt_seq}_b38
–tss_distance: distance between variant and transcription start site (TSS). Positive when variant is downstream of the TSS, negative otherwise
–chr: chromosome (for the variant, same as gene_chr for cis-eQTLs)
–variant_pos: position of the first reference base of the variant
–ref: reference sequence of the variant
–alt: alternate sequence of the variant
–num_alt_per_site: number of alternative alleles observed at this site
–rs_id_dbSNP151_GRCh38p7: dbSNP151 rsID
–minor_allele_samples: number of samples carrying the minor allele
–minor_allele_count: total number of minor alleles across individuals
–maf: minor allele frequency observed in the set of donors for a given tissue
–ref_factor: ‘1’, when the minor allele is the alt base, ‘-1’ when the minor allele is the reference base
–pval_nominal: nominal p-value associated with the most significant variant for this gene
–slope: regression slope
–slope_se: standard error of the regression slope
–pval_perm: permutation p-value
–pval_beta: beta-approximated permutation p-value
–qval: Storey q-value derived from pval_beta
–pval_nominal_threshold: nominal p-value threshold for calling a variant-gene pair significant for the gene
–log2_aFC: allelic fold change (aFC) in log2 scale (see Mohammadi et al., Genome Res. 2017)
–log2_aFC_lower: lower bound of the 95% confidence interval of log2(aFC)
–log2_aFC_upper: upper bound of the 95% confidence interval of log2(aFC)
–Significant variant-gene pairs
–File extensions: .signif_variant_gene_pairs.txt.gz, .sqtl_signifpairs.txt.gz
–Column headers:
–variant_id: variant ID in the format {chr}{pos_first_ref_base}{ref_seq}_{alt_seq}_b38
–gene_id: GENCODE/Ensembl gene ID
–tss_distance: distance between variant and transcription start site. Positive when variant is downstream of the TSS, negative otherwise
–ma_samples: number of samples carrying the minor allele
–ma_count: total number of minor alleles across individuals
–maf: minor allele frequency observed in the set of donors for a given tissue
–pval_nominal: nominal p-value
–slope: regression slope
–slope_se: standard error of the regression slope
–pval_nominal_threshold: nominal p-value threshold for calling a variant-gene pair significant for the gene
–min_pval_nominal: smallest nominal p-value for the gene
–pval_beta: beta-approximated permutation p-value for the gene
–All variant-gene pairs tested
–File extension: .allpairs.txt.gz, .sqtl_allpairs.txt.gz
–Column headers:
–gene_id: GENCODE/Ensembl gene ID
–variant_id: variant ID in the format {chr}{pos_first_ref_base}{ref_seq}_{alt_seq}_b38
–tss_distance: distance between variant and transcription start site. Positive when variant is downstream of the TSS, negative otherwise
–ma_samples: number of samples carrying the minor allele
–ma_count: total number of minor alleles across individuals
–pval_nominal: nominal p-value
–slope: regression slope
–slope_se: standard error of the regression slope
# GTEx EUR pairs in Genome Reference Consortium Human Build 38 (GRCh38)
# Position from GRCh37 needed to harmonize with the ASD and 1000 Genomes data in GRCh37
# https://www.nature.com/articles/s41591-021-01310-z
# Get the GTEx build-converting file
<!-- wget https://storage.googleapis.com/gtex_analysis_v8/reference/GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt.gz -->
<!-- mv GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt.gz /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR -->
<!-- gunzip GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt.gz -->
<!-- # Get GTEx files -->
<!-- wget https://storage.googleapis.com/gtex_analysis_v8/single_tissue_qtl_data/GTEx_Analysis_v8_eQTL_EUR.tar -->
<!-- gunzip GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt.gz -->
<!-- # Uncompress them: they include egene, pairs, and other files -->
<!-- tar -xvf GTEx_Analysis_v8_eQTL_EUR.tar -->
<!-- # Create "egene" and "pairs" directories -->
<!-- # mkdir egene -->
<!-- cd eqtls -->
<!-- cp *egenes.txt.gz /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene -->
<!-- # Unzip the tar files once they are in the right directories -->
<!-- gunzip *.gz -->
<!-- mkdir pairs -->
<!-- cd eqtls -->
<!-- cp *signif_pairs.txt.gz /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs -->
<!-- # Unzip the tar files once they are in the right directories -->
<!-- gunzip *.gz -->–Qi et al. brain eQTL summary data
–GTEx-brain eQTL data (estimated effective n = 233)
–GTEx-brain eQTL summary data (Qi et al. 2018 Nat Commun) in SMR binary (BESD) format:
–GTEx-brain.tar.gz (1.1 GB)
–This is a set of eQTL data from a meta-analysis of 10 brain regions in GTEx v6
–(GTEx Consortium 2017 Nature) correcting for sample overlap by the MeCS method.
–Only SNPs within 1Mb distance from each probe are available.
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/brain -->
<!-- tar xvfz GTEx-brain.tar.gz -->
<!-- # Copy over the ASD and 1000 Genomes binary files -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala -->
<!-- cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/brain -->
<!-- cp GTEx_Brain_std.besd GTEx_Brain_std.epi GTEx_Brain_std.esi /n/holylfs/LABS/lemos_lab/Users/cdadams/brain -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/brain -->
–Brain-eMeta eQTL data (estimated effective n = 1194)
–Brain-eMeta eQTL summary data (Qi et al. 2018 Nat Commun) in SMR binary (BESD) format:
–Brain-eMeta.tar.gz (1.1 GB)
–This is a set of eQTL data from a meta-analysis of GTEx brain (GTEx Consortium 2017 Nature), CMC (Fromer et al. 2016 Nat Neurosci), and ROSMAP (Ng et al. 2017 Nat Neurosci) by MeCS. Only SNPs within 1Mb distance from each probe are available.
<!-- wget https://cnsgenomics.com/data/SMR/Brain-eMeta.tar.gz -->
<!-- tar xvfz Brain-eMeta.tar.gz -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/Brain_eMeta/Brain-eMeta -->
<!-- cp Brain-eMeta.besd Brain-eMeta.epi Brain-eMeta.esi /n/holylfs/LABS/lemos_lab/Users/cdadams/Brain_eMeta/ -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/brain/ -->
<!-- cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/Brain_eMeta -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/Brain_eMeta/ -->
–The 1000 Genomes LD (bed, bim, fam) files are used by SMR
–Get the European LD reference panel binaries from 1000 genomes
# devtools::install_github("explodecomputer/genetics.binaRies")
# genetics.binaRies::get_plink_binary()
# https://mrcieu.github.io/ieugwasr/articles/local_ld.html<!-- wget http://fileserve.mrcieu.ac.uk/ld/1kg.v3.tgz -->
<!-- tar -zxvf 1kg.v3.tgz -->
<!-- # Copy the LD reference files over: they need to be in every directory SMR is run from... -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/smr_eqtlGen -->
<!-- cp EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR -->
–Need to know which tissue the files are for to merge them into one
–Add a column “name” to end of each file and populate it with the filename
–https://stackoverflow.com/questions/68039474/add-column-of-filename-to-hundreds-of-files-bash
<!-- tmp=$(mktemp) || { ret="$?"; printf 'Failed to create temp file\n'; exit "$ret"; } -->
<!-- for file in *.txt; do -->
<!-- awk 'BEGIN{OFS="\t"} {print $0, (FNR>1 ? FILENAME : "name")}' "$file" > "$tmp" && -->
<!-- mv -- "$tmp" "$file" || exit -->
<!-- done -->
# The files are big; otherwise I'd have just read the files into R and merged them
# Merge all into a single file keeping only first header
# https://unix.stackexchange.com/questions/60577/concatenate-multiple-files-with-same-header
<!-- awk ' -->
<!-- FNR==1 && NR!=1 { while (/^<header>/) getline; } -->
<!-- 1 {print} -->
<!-- ' *.txt >all_pairs.txt --><!-- ssh cdadams@login.rc.fas.harvard.edu -->
<!-- salloc -n 1 -p shared -t 0-09:00 --mem=9000 -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR -->
#library(bannerCommenter)
library(biomaRt)
library(Cairo)
library(car)
#library(caret)
#library(CCA)
#library(CCP)
#library(corrplot)
#library(CMplot)
#library(dagitty)
library(data.table)
library(DataCombine)
library(dplyr)
library(DT)
library(e1071)
#library(EnhancedVolcano)
library(EnsDb.Hsapiens.v86)
# library(ensembldb)
# library(foreach)
# library(GGally)
#library(GeneCorrAtlas)
#library(ggdag)
library(ggforestplot)
library(ggplot2)
library(ggpubr)
library(ggrepel)
# library(ggVennDiagram)
# library(glmnet)
library(gplots)
library(graphics)
library(grid)
library(gridExtra)
library(Hmisc)
#library(htmltools)
#library(hyprcoloc)
#library(igraph)
library(kableExtra)
library(lme4)
#library(magick)
library(MASS)
library(meta)
#library(metaviz)
#library(MRInstruments)
#library(oncofunco)
#library(packcircles)
#library(parallel)
library(pheatmap)
library(plyr)
library(qqman)
library(RadialMR)
library(RColorBrewer)
library(rcompanion)
library(reshape2)
#library(simex)
library(sqldf)
library(strex)
library(stringr)
library(tidyr)
library(tidyverse)
library(tm)
library(TwoSampleMR)
library(vcd)
library(viridis)
#install.packages("glmnet", repos = "https://cran.us.r-project.org")
# Cleaned GTEx file used across all analyses:
# gtex_pairs_rsid_thresh <- readRDS('/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid_thresh.Rda')# Table for converting from b38 to b37
# Obtained from: https://www.nature.com/articles/s41591-021-01310-z
tbl=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/GTEx_Analysis_2017-06-05_v8_WholeGenomeSeq_838Indiv_Analysis_Freeze.lookup_table.txt")
tbl=as.data.frame(tbl)
# Examine an RP SNP to see if the position is different for GRCh38 vs GRCh37: it is.
# rsids are stable but position isn't between builds.
# Need to use the GRCh37 position for the autism data, as the LD reference panel
# and the autism GWAS are GRCh37
# Regular MR doesn't use position but SMR does, which is why I didn't have a problem
# when running earlier versions of this in 2020 for the RP meta-analyses in brain.
# All eqtls
setwd("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs")
gtex_pairs=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/all_pairs.txt")#pairs
# Get the rsids for gtex_pairs
gtex_pairs_rsid=merge(gtex_pairs, tbl, by="variant_id", all.x=TRUE)
# Get the gene names
ensemblsIDS=gtex_pairs_rsid$phenotype_id
# Remove the version numbers at the end of the ensembl names
gtex_pairs_rsid$GENEID <- gsub('\\..+$', '', gtex_pairs_rsid$phenotype_id)
# Fetch the gene names: convert from ensembl.gene to gene.symbol
ensembl.genes=as.character(gtex_pairs_rsid$GENEID)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype =
"GENEID", columns = c("SYMBOL","GENEID"))
# Merge in the gene names
gtex_pairs_rsid2 <- merge(x = geneIDs1,
y = gtex_pairs_rsid,
by.x="GENEID",
by.y="GENEID")
# Get the names for the tissues nice (remove .txt and hyphens)
gtex_pairs_rsid2$tissue=sub("\\..*", "", gtex_pairs_rsid2$name)
gtex_pairs_rsid2$tissue2=gsub("_", " ", gtex_pairs_rsid2$tissue)
# Convert to numeric
cols.num <- c("tss_distance","ma_samples",
"ma_count","maf",
"pval_nominal","slope","slope_se",
"min_pval_nominal","pval_beta",
"pval_nominal_threshold")
gtex_pairs_rsid2[cols.num] <- sapply(gtex_pairs_rsid2[cols.num],as.numeric)
# Calculate z-scores for SMR
gtex_pairs_rsid2$z=qnorm(gtex_pairs_rsid2$pval_nominal)
# Calculate adjusted SE for SMR using z-score,
# per: https://cnsgenomics.com/software/smr/#BESDformat
# SE = b / z*
# However, SE isn't read in with GTEx fastqtl-nominal-format
gtex_pairs_rsid2$smr_SE=abs(gtex_pairs_rsid2$slope/gtex_pairs_rsid2$z)
# Extract the CHR37 position
# Do in two steps with 'sub'
# Extract everything after first hyphen
gtex_pairs_rsid2$position_37a=sub(".*?_",'', gtex_pairs_rsid2$variant_id_b37)
# Extract everything before first hyphen: didn't work!
#gtex_pairs_rsid2$position_37=sub("\\_.*", "", gtex_pairs_rsid2$position_37a)
#gtex_pairs_rsid2$chr_num=substr(gtex_pairs_rsid2$chr, 4, 4)
#head(gtex_pairs_rsid2)
# Get the P<5e-04 data (before we used P<5e-06, but others use P<5e-04)
gtex_pairs_rsid_thresh=gtex_pairs_rsid2[which(gtex_pairs_rsid2$pval_nominal<5e-04),]
# Get chromosome number
gtex_pairs_rsid_thresh$chr_num=str_remove(gtex_pairs_rsid_thresh$chr , "chr")
# Look at position_37
#is.na(gtex_pairs_rsid_thresh$position_37)
gtex_pairs_rsid_thresh_complete=gtex_pairs_rsid_thresh[complete.cases(gtex_pairs_rsid_thresh[ , 28]),]
gtex_pairs_rsid_thresh=gtex_pairs_rsid_thresh_complete
# Save gtex_pairs_rsid_thresh and gtex_pairs_rsid2 as .Rda files individually
saveRDS(gtex_pairs_rsid_thresh, file = '/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid_thresh.Rda', compress = FALSE)
#saveRDS(gtex_pairs_rsid2, file = '/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid2.Rda', compress = FALSE)<!-- setwd("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs") -->
<!-- table(gtex_pairs_rsid_thresh$tissue) -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/ -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Visceral_Omentum -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adrenal_Gland -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Artery_Aorta -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Artery_Coronary -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Artery_Tibial -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Anterior_cingulate_cortex_BA24 -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Caudate_basal_ganglia -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellar_Hemisphere -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellum -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cortex -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Frontal_Cortex_BA9 -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hippocampus -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hypothalamus -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Nucleus_accumbens_basal_ganglia -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Putamen_basal_ganglia -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Spinal_cord_cervical_c1 -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Substantia_nigra -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Breast_Mammary_Tissue -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Cells_Cultured_fibroblasts -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Cells_EBV_transformed_lymphocytes -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Colon_Sigmoid -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Colon_Transverse -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Esophagus_Gastroesophageal_Junction -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Esophagus_Mucosa -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Esophagus_Muscularis -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Heart_Atrial_Appendage -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Heart_Left_Ventricle -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Kidney_Cortex -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Liver -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Lung -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Minor_Salivary_Gland -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Muscle_Skeletal -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Nerve_Tibial -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Ovary -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Pituitary -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Pancreas -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Prostate -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Skin_Not_Sun_Exposed_Suprapubic -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Skin_Sun_Exposed_Lower_leg -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Spleen -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Stomach -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Testis -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Uterus -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Vagina -->
<!-- # mkdir /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Whole_Blood -->
gtex_pairs_rsid_thresh <- readRDS('/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid_thresh.Rda')
split <- split(gtex_pairs_rsid_thresh, gtex_pairs_rsid_thresh$tissue)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
files <- mget(ls())
length(files)
regions=c("Brain_Amygdala","Brain_Anterior_cingulate_cortex_BA24",
"Brain_Caudate_basal_ganglia","Brain_Cerebellar_Hemisphere",
"Brain_Cerebellum","Brain_Cortex",
"Brain_Frontal_Cortex_BA9","Brain_Hippocampus",
"Brain_Hypothalamus", "Brain_Nucleus_accumbens_basal_ganglia",
"Brain_Putamen_basal_ganglia","Brain_Spinal_cord_cervical_c-1",
"Brain_Substantia_nigra")
rm(gtex_pairs_rsid_thresh)
rm(`Cells_EBV-transformed_lymphocytes`)
rm(`Brain_Spinal_cord_cervical_c-1`)
rm(split)
rm(i)
rm(regions)
# Had to change the names for tissues with weird symbols
Cells_EBV_transformed_lymphocytes=`Cells_EBV-transformed_lymphocytes`
Brain_Spinal_cord_cervical_c1=`Brain_Spinal_cord_cervical_c-1`
# Save .csvs
write.csv(Brain_Amygdala,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala/Brain_Amygdala.csv")
write.csv(Brain_Anterior_cingulate_cortex_BA24,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Anterior_cingulate_cortex_BA24/Brain_Anterior_cingulate_cortex_BA24.csv")
write.csv(Brain_Caudate_basal_ganglia,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Caudate_basal_ganglia/Brain_Caudate_basal_ganglia.csv")
write.csv(Brain_Cerebellar_Hemisphere,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellar_Hemisphere/Brain_Cerebellar_Hemisphere.csv")
write.csv(Brain_Cerebellum,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellum/Brain_Cerebellum.csv")
write.csv(Brain_Cortex,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cortex/Brain_Cortex.csv")
write.csv(Brain_Frontal_Cortex_BA9,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Frontal_Cortex_BA9/Brain_Frontal_Cortex_BA9.csv")
write.csv(Brain_Hippocampus,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hippocampus/Brain_Hippocampus.csv")
write.csv(Brain_Hypothalamus,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hypothalamus/Brain_Hypothalamus.csv")
write.csv(Brain_Nucleus_accumbens_basal_ganglia,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Nucleus_accumbens_basal_ganglia/Brain_Nucleus_accumbens_basal_ganglia.csv")
write.csv(Brain_Putamen_basal_ganglia,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Putamen_basal_ganglia/Brain_Putamen_basal_ganglia.csv")
write.csv(Brain_Spinal_cord_cervical_c1,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Spinal_cord_cervical_c1/Brain_Spinal_cord_cervical_c1.csv")
write.csv(Brain_Substantia_nigra,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Substantia_nigra/Brain_Substantia_nigra.csv")
write.csv(Vagina,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Vagina/Vagina.csv")
write.csv(Adipose_Subcutaneous,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Adipose_Subcutaneous.csv")
write.csv(Adipose_Visceral_Omentum,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Adipose_Visceral_Omentum.csv")
write.csv(Adrenal_Gland,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Adrenal_Gland.csv")
write.csv(Artery_Aorta,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Artery_Aorta.csv")
write.csv(Artery_Coronary,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Artery_Coronary.csv")
write.csv(Artery_Tibial,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Artery_Tibial.csv")
write.csv(Breast_Mammary_Tissue,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Breast_Mammary_Tissue.csv")
write.csv(Cells_Cultured_fibroblasts,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Cells_Cultured_fibroblasts.csv")
write.csv(Cells_EBV_transformed_lymphocytes,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Cells_EBV_transformed_lymphocytes.csv")
write.csv(Colon_Sigmoid,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Colon_Sigmoid.csv")
write.csv(Colon_Transverse,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Colon_Transverse.csv")
write.csv(Esophagus_Gastroesophageal_Junction,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Esophagus_Gastroesophageal_Junction.csv")
write.csv(Esophagus_Mucosa,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Esophagus_Mucosa.csv")
write.csv(Esophagus_Muscularis,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Esophagus_Muscularis.csv")
write.csv(Heart_Atrial_Appendage,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Heart_Atrial_Appendage.csv")
write.csv(Heart_Left_Ventricle,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Heart_Left_Ventricle.csv")
write.csv(Kidney_Cortex,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Kidney_Cortex.csv")
write.csv(Liver,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Liver.csv")
write.csv(Lung,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Lung.csv")
write.csv(Minor_Salivary_Gland,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Minor_Salivary_Gland.csv")
write.csv(Muscle_Skeletal,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Muscle_Skeletal.csv")
write.csv(Nerve_Tibial,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Nerve_Tibial.csv")
write.csv(Ovary,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Ovary.csv")
write.csv(Pancreas,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Pancreas.csv")
write.csv(Prostate,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Prostate.csv")
write.csv(Skin_Not_Sun_Exposed_Suprapubic,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Skin_Not_Sun_Exposed_Suprapubic.csv")
write.csv(Skin_Sun_Exposed_Lower_leg,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Skin_Sun_Exposed_Lower_leg.csv")
write.csv(Spleen,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Spleen.csv")
write.csv(Stomach,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Stomach.csv")
write.csv(Testis,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Testis.csv")
write.csv(Thyroid,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Thyroid.csv")
write.csv(Uterus,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Uterus.csv")
write.csv(Vagina,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Vagina.csv")
write.csv(Whole_Blood,"/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous/Whole_Blood.csv")
# Copy the LD binaries and ASD data into the respective tissue folders for the pairs
# cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/egene/Brain_Amygdala
# cp smr_asd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala/
# cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala/
#
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Amygdala
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Anterior_cingulate_cortex_BA24
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Caudate_basal_ganglia
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellar_Hemisphere
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cerebellum
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Cortex
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Frontal_Cortex_BA9
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hippocampus
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Hypothalamus
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Nucleus_accumbens_basal_ganglia
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Putamen_basal_ganglia
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Spinal_cord_cervical_c1
# cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Brain_Substantia_nigra<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR -->
<!-- mkdir all_tissue_SMR -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/Adipose_Subcutaneous -->
<!-- cp smr_asd EUR.bed EUR.bim EUR.fam /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR -->
<!-- cp /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid_thresh.Rda /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR -->
<!-- cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR -->
<!-- ls -->#cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR
gtex_pairs_rsid_thresh <- readRDS('/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/pairs/gtex_pairs_rsid_thresh.Rda')
gtex_pairs_rsid_thresh$symbol_tissue=paste(gtex_pairs_rsid_thresh$SYMBOL, gtex_pairs_rsid_thresh$tissue)
gtex_pairs_rsid_thresh$symbol_tissue <- gsub(" ", "_", gtex_pairs_rsid_thresh$symbol_tissue)
# gtex_pairs_rsid_thresh$symbol_tissue_rsid=paste(gtex_pairs_rsid_thresh$SYMBOL, gtex_pairs_rsid_thresh$tissue, gtex_pairs_rsid_thresh$rs_id_dbSNP151_GRCh38p7)
# gtex_pairs_rsid_thresh$symbol_tissue_rsid <- gsub(" ", "_", gtex_pairs_rsid_thresh$symbol_tissue_rsid)
all_smr=gtex_pairs_rsid_thresh
rm(gtex_pairs_rsid_thresh)
myvars=c("symbol_tissue","rs_id_dbSNP151_GRCh38p7","tss_distance", "pval_nominal", "slope",
"SYMBOL", "alt", "ref", "position_37", "chr_num")
all_smr=gtex_pairs_rsid_thresh[myvars]
# Convert dot in rsid to NA and remove NAs
all_smr$rs_id_dbSNP151_GRCh38p7[all_smr$rs_id_dbSNP151_GRCh38p7 == "."] <- NA
all_smr <- na.omit(all_smr)
# Remove rows with more than one rsid
#complete %>% "["(.,14783:14783,)
all_smr=all_smr[- grep(",", all_smr$rs_id_dbSNP151_GRCh38p7),]
all_smr=all_smr[order(all_smr$rs_id_dbSNP151_GRCh38p7),]
# For making besd file: use "symbol_tissue" for SMR and "symbol_tissue_rsid" for MR
myvars=c("symbol_tissue","rs_id_dbSNP151_GRCh38p7","tss_distance", "pval_nominal", "slope")
#myvars=c("symbol_tissue_rsid","rs_id_dbSNP151_GRCh38p7","tss_distance", "pval_nominal", "slope")
all_smr_besd=all_smr[myvars]
all_smr_besd=all_smr_besd[order(all_smr_besd$rs_id_dbSNP151_GRCh38p7),]
all_smr_besd$id <- 1:nrow(all_smr_besd)
all_smr=all_smr[order(all_smr$rs_id_dbSNP151_GRCh38p7),]
all_smr$id <- 1:nrow(all_smr)
head(all_smr)
all_smr_besd$id=NULL
# Deal with duplicates -----------------------------------------
dup=all_smr_besd[duplicated(all_smr_besd[,1:2]),]
all_smr_besd=all_smr_besd[!duplicated(all_smr_besd[c(1,2)]),]
# Save without column names for making the besd, epi, and esi files for *SMR*: pay attention!
write.table(all_smr_besd,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/all_smr.txt",
sep = "\t",
row.names = FALSE, col.names = FALSE, quote=FALSE)
saveRDS(all_smr, file = '/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/all_smr_FULL_for_epi.Rda', compress = FALSE)# Remove the quotes
<!-- sed 's/\"//g' all_smr.txt > all_smr -->
<!-- # SMR command for making the besd, epi, and esi files: *SMR* -->
<!-- smr_Linux smr --eqtl-summary all_smr --fastqtl-nominal-format --make-besd --out mybesd_all_smr -->
# SMR command for making the besd, epi, and esi files: *MR*: moved to all_smr_res_MR_pay_attention.sh and
# script_for_besd_MR_not_SMR_ASD.R
#smr_Linux smr --eqtl-summary all_smr_MR --fastqtl-nominal-format --make-besd --out mybesd_all_smr_MR
epi=read.table("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/mybesd_all_smr.epi")
# Columns are chromosome, probe ID(can be the ID of an exon or a transcript for RNA-seq data),
# genetic distance (can be any arbitary value), physical position, gene ID and gene orientation
epi$id <- 1:nrow(epi)
colnames(epi)[1]="chr_num"
colnames(epi)[2]="symbol_tissue"
colnames(epi)[3]="genetic distance"
colnames(epi)[4]="position"
colnames(epi)[5]="SYMBOL2"
colnames(epi)[6]="gene orientation"
outepi=merge(epi, all_smr, by="symbol_tissue")
outepi=outepi[order(outepi$id.x),]
z <- outepi[!duplicated(outepi$symbol_tissue),]
z=z[order(z$id.x), ]
z$symbol_tissue2=z$symbol_tissue
z$positionDummy=z$symbol_tissue
# Columns are:
# chromosome,
# probe ID (can be the ID of an exon or a transcript for RNA-seq data),
# genetic distance (can be any arbitary value),
# physical position,
# gene ID
# and gene orientation
myvars=c("chr_num.y","symbol_tissue","genetic distance","position_37","SYMBOL","id.x")
all_smr_epi=z[myvars]
# Save for *SMR*
write.table(all_smr_epi,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/mybesd_all_smr.epi",
sep = "\t",
row.names = FALSE, col.names = FALSE, quote=FALSE)# The esi file needs to be hand-modified to include chr, rsid, genetic distance (dummy),
# position, alt, ref, EAF
# The order of the rsids needs to be maintained, as the headers are removed and they must
# match the order in the besd file
esi=read.table("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/mybesd_all_smr.esi")
colnames(esi)[1]="chr_num"
colnames(esi)[2]="rs_id_dbSNP151_GRCh38p7"
colnames(esi)[3]="genetic distance"
colnames(esi)[4]="position"
colnames(esi)[5]="alt"
colnames(esi)[6]="ref"
colnames(esi)[7]="EAF"
# Must maintain the order, as it matches with besd
esi$id <- 1:nrow(esi)
out=merge(esi, all_smr, by="rs_id_dbSNP151_GRCh38p7")
out=out[order(out$id.x), ]
myvars=c("chr_num.y","rs_id_dbSNP151_GRCh38p7","genetic distance","position_37","alt.y","ref.y","EAF", "id.x")
all_smr_esi=out[myvars]
clean=all_smr_esi %>% distinct(rs_id_dbSNP151_GRCh38p7, .keep_all = TRUE)
clean=clean[order(clean$id.x), ]
clean$id.x=NULL
# Save for SMR
write.table(clean,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/mybesd_all_smr.esi",
sep = "\t",
row.names = FALSE, col.names = FALSE, quote=FALSE)# Ran this with a bash script: important to have the "out" (.o) file to know if it completed
#smr_Linux smr --bfile EUR --gwas-summary smr_asd --beqtl-summary mybesd_all_smr --peqtl-smr 5e-4 --out mysmr_all_smr --thread-num 10
# PAY ATTENTION: save all results: /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/
smr_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/mysmr_all_smr.smr")
# rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
smr_asd$tissue=sub('.-*', '-', smr_asd$probeID)
smr_asd$tissue=str_after_nth(smr_asd$probeID, "_", 1)
smr_asd=as.data.frame(smr_asd)
smr_asd=as.data.frame(smr_asd)
dim(smr_asd)
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_smr_asd=smr_asd[which(smr_asd$Gene %in% combo4$gene),]
rp_smr_asd=rp_smr_asd[order(rp_smr_asd$p_SMR),]
rp_smr_asd=as.data.frame(rp_smr_asd)
write.table(rp_smr_asd,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/rp_asd.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
rp_smr_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/rp_asd.txt")
rp_smr_asd=as.data.frame(rp_smr_asd)
# Read-in eMeta results
rp_emeta_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/rp_mysmr_Brain_Meta-analysis_(eMeta).txt")
rp_emeta_asd=as.data.frame(rp_emeta_asd)
rp_emeta_asd$tissue="eMeta_(cross-brain)"
rp_emeta_asd$fdr=NULL
rp_emeta_asd=rp_emeta_asd[order(rp_emeta_asd$p_SMR),]
# Read-in GTEx brain results
rp_gtex_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/rp_mysmr_Brain_Meta-analysis_(GTEx).txt")
rp_gtex_asd$fdr=NULL
rp_gtex_asd$tissue="GTEx_(cross-brain)"
rp_gtex_asd=as.data.frame(rp_gtex_asd)
rp_gtex_asd=rp_gtex_asd[order(rp_gtex_asd$p_SMR),]
# Read-in eqtlgen results
rp_eqtlgen_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/rp_mysmr_eqtlGen_Whole_Blood_test.txt")
rp_eqtlgen_asd$tissue="Whole_Blood_(eqtlGen)"
rp_eqtlgen_asd$fdr=NULL
rp_eqtlgen_asd$GENEID=NULL
rp_eqtlgen_asd$SYMBOL=NULL
# rbind GTEx tissues, eMeta, GTEx brain, and eqtlGen results
WB_asd=rbind(rp_smr_asd,rp_emeta_asd,rp_gtex_asd,rp_eqtlgen_asd)
WB_asd$GWAS="ASD"
write.table(WB_asd,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/RP_WB_asd.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)–Alzheimer’s disease
–Dataset: ieu-a-297
–PMID 24162737
–Year 2013
–Category Disease
–Sub category Psychiatric / neurological
–Population European
–Sex Males and Females
–ncase 17,008
–ncontrol 37,154
–Sample size 54,162
–Number of SNPs 7,055,882
–Unit log odds
–Priority 2
–Author Lambert
–Consortium IGAP
–Ontology NA
–Build HG19/GRCh37
–Note Effect allele frequencies are missing; forward(+) strand
<!-- wget https://gwas.mrcieu.ac.uk/files/ieu-a-297/ieu-a-297.vcf.gz --># Read two times the vcf file, first for the columns names, second for the data
tmp_vcf<-readLines(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/ieu-a-297.vcf.gz"))
tmp_vcf_data<-read.table(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/ieu-a-297.vcf.gz"), stringsAsFactors = FALSE)
# Filter for the columns names
tmp_vcf<-tmp_vcf[-(grep("#CHROM",tmp_vcf)+1):-(length(tmp_vcf))]
vcf_names<-unlist(strsplit(tmp_vcf[length(tmp_vcf)],"\t"))
names(tmp_vcf_data)<-vcf_names
head(tmp_vcf_data) #ES:SE:LP:AF:SS:ID
colnames(tmp_vcf_data)[1]="chromosome"
colnames(tmp_vcf_data)[2]="POS"
colnames(tmp_vcf_data)[3]="SNP"
colnames(tmp_vcf_data)[4]="A2"
colnames(tmp_vcf_data)[5]="A1"
colnames(tmp_vcf_data)[10]="mix"
tmp_vcf_data$b=sub("\\:.*", "", tmp_vcf_data$mix)
tmp_vcf_data$SE.1=str_before_nth(tmp_vcf_data$mix, ":", 2)
tmp_vcf_data$se=sub('.*:', '', tmp_vcf_data$SE.1)
tmp_vcf_data$p.1=str_before_nth(tmp_vcf_data$mix, ":", 3)
tmp_vcf_data$p=sub('.*:', '', tmp_vcf_data$p.1)
# tmp_vcf_data$eaf.1=str_before_nth(tmp_vcf_data$mix, ":", 4)
# tmp_vcf_data$eaf=sub('.*:', '', tmp_vcf_data$eaf.1)
# tmp_vcf_data$N.1=str_before_nth(tmp_vcf_data$mix, ":", 5)
# tmp_vcf_data$N=sub('.*:', '', tmp_vcf_data$N.1)
tmp_vcf_data$p=as.numeric(as.character(tmp_vcf_data$p))
tmp_vcf_data$b=as.numeric(as.character(tmp_vcf_data$b))
#tmp_vcf_data$N=as.numeric(as.character(tmp_vcf_data$N))
#tmp_vcf_data$eaf=as.numeric(as.character(tmp_vcf_data$eaf))
tmp_vcf_data$se=as.numeric(as.character(tmp_vcf_data$se))
summary(tmp_vcf_data$p)
tmp_vcf_data$pval=10^(-tmp_vcf_data$p)
tmp_vcf_data$pval2=10^(tmp_vcf_data$p)
summary(tmp_vcf_data$pval)
summary(tmp_vcf_data$p)
tmp_vcf_data$p=tmp_vcf_data$pval
tmp_vcf_data$n=54162
tmp_vcf_data$freq=NA
# Select asd variables for SMR
myvars=c("SNP", "A1", "A2", "freq", "b", "se", "p", "n" )
smr_alzheimers=tmp_vcf_data[myvars]
smr_alzheimers=smr_alzheimers[order(smr_alzheimers$p),]
alz_sigs=smr_alzheimers[which(smr_alzheimers$p<0.00000005),]
write.table(smr_alzheimers, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/smr_alzheimers.txt",
sep="\t", row.names = FALSE, col.names = TRUE)# Remove quotes to run in SMR
<!-- sed 's/\"//g' smr_alzheimers.txt > smr_alzheimers -->
# SMR of eqtlGen, GTEx brain, and eMeta brain on Alzheimers
# Ran with a bash (slurm / sbatch) script
# smr_Linux smr --bfile EUR --gwas-summary smr_alzheimers --beqtl-summary cis-eQTLs-full_eQTLGen_AF_incl_nr_formatted_20191212.new.txt_besd-dense --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_eqtlgen_alzheimers --thread-num 10
# smr_Linux smr --bfile EUR --gwas-summary smr_alzheimers --beqtl-summary Brain-eMeta --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_Brain_eMeta_alzheimers --thread-num 10
# smr_Linux smr --bfile EUR --gwas-summary smr_alzheimers --beqtl-summary GTEx_Brain_std --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_Brain_GTEx_alzheimers --thread-num 10
setwd("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers")
alz_smr=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_all_smr_alzheimers_SR_alzOnly.smr")
# rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
alz_smr$tissue=sub('.-*', '-', alz_smr$probeID)
alz_smr$tissue=str_after_nth(alz_smr$probeID, "_", 1)
alz_smr=as.data.frame(alz_smr)
head(alz_smr)# gene.alz=table(alz_smr$Gene)
# dim(gene.alz)
# bon.alz=0.05/31293
# alz.bon=alz_smr[which(alz_smr$p_SMR<bon.alz),]
# dim(alz.bon)
# alz_smr=alz_smr[order(alz_smr$p_SMR),]
# head(alz_smr$Gene, n=20)
# head(alz_smr)
# apob=alz_smr[which(alz_smr$Gene=="APOB"),]# Select the NUC set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_alz_smr=alz_smr[which(alz_smr$Gene %in% combo4$gene),]
rp_alz_smr=rp_alz_smr[order(rp_alz_smr$p_SMR),]
write.table(rp_alz_smr, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/rp_alz_smr.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)# Read-in emeta brain results
alz_emeta=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_Brain_eMeta_alzheimers.smr")
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_alz_emeta=alz_emeta[which(alz_emeta$Gene %in% combo4$gene),]
rp_alz_emeta=rp_alz_emeta[order(rp_alz_emeta$p_SMR),]
rp_alz_emeta=as.data.frame(rp_alz_emeta)
rp_alz_emeta$tissue="eMeta_(cross-brain)"
write.table(rp_alz_emeta, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/rp_alz_emeta.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
# Read-in GTEx brain results
alz_gtex=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_Brain_GTEx_alzheimers.smr")
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_alz_gtex=alz_gtex[which(alz_gtex$Gene %in% combo4$gene),]
rp_alz_gtex=rp_alz_gtex[order(rp_alz_gtex$p_SMR),]
rp_alz_gtex=as.data.frame(rp_alz_gtex)
rp_alz_gtex$tissue="GTEx_(cross-brain)"
write.table(rp_alz_gtex, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/rp_alz_gtex.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
# Read-in eqtlgen results
alz_eqtlgen=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/mysmr_eqtlgen_alzheimers.smr")
# Get the gene symbols
ensemblsIDS=alz_eqtlgen$Gene
alz_eqtlgen$GENEID=alz_eqtlgen$Gene
# Remove the version numbers at the end of the ensembl names
#gtex_pairs_rsid$GENEID <- gsub('\\..+$', '', gtex_pairs_rsid$phenotype_id)
# Fetch the gene names: convert from ensembl.gene to gene.symbol
ensembl.genes=as.character(alz_eqtlgen$GENEID)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype =
"GENEID", columns = c("SYMBOL","GENEID"))
# Merge in the gene names
alz_eqtlgen2 <- merge(x = geneIDs1,
y = alz_eqtlgen,
by.x="GENEID",
by.y="GENEID")
alz_eqtlgen2$GENEID=NULL
alz_eqtlgen2$Gene=alz_eqtlgen2$SYMBOL
alz_eqtlgen2$SYMBOL=NULL
alz_eqtlgen=alz_eqtlgen2
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_alz_eqtlgen=alz_eqtlgen[which(alz_eqtlgen$Gene %in% combo4$gene),]
rp_alz_eqtlgen=rp_alz_eqtlgen[order(rp_alz_eqtlgen$p_SMR),]
rp_alz_eqtlgen=as.data.frame(rp_alz_eqtlgen)
rp_alz_eqtlgen$tissue="Whole_Blood_(eqtlGen)"
write.table(rp_alz_eqtlgen, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/rp_alz_eqtlgen.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
WB_alz=rbind(rp_alz_smr,rp_alz_emeta,rp_alz_gtex,rp_alz_eqtlgen)
WB_alz$GWAS="Alzheimers"
write.table(WB_alz, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/RP_WB_alz.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
WB_alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/RP_WB_alz.txt")
write.table(WB_alz,
file = "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/rp_alzheimers.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)–Parkinson’s disease
–Dataset: ieu-b-7
–Year 2019
–Category Binary
–Sub category Psychiatric / neurological
–Population European
–Sex Males and Females
–ncase 33,674
–ncontrol 449,056
–Sample size 482,730
–Number of SNPs 17,891,936
–Unit NA
–Priority 1
-Author Nalls MA
–Consortium International Parkinson’s Disease Genomics Consortium
–Ontology NA
–Build HG19/GRCh37
<!-- wget https://gwas.mrcieu.ac.uk/files/ieu-b-7/ieu-b-7.vcf.gz --># Read two times the vcf file, first for the columns names, second for the data
tmp_vcf<-readLines(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/ieu-b-7.vcf.gz"))
tmp_vcf_data<-read.table(gzfile("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/ieu-b-7.vcf.gz"), stringsAsFactors = FALSE)
# filter for the columns names
tmp_vcf<-tmp_vcf[-(grep("#CHROM",tmp_vcf)+1):-(length(tmp_vcf))]
vcf_names<-unlist(strsplit(tmp_vcf[length(tmp_vcf)],"\t"))
names(tmp_vcf_data)<-vcf_names
head(tmp_vcf_data) #ES:SE:LP:AF:SS:ID
colnames(tmp_vcf_data)[1]="chromosome"
colnames(tmp_vcf_data)[2]="POS"
colnames(tmp_vcf_data)[3]="SNP"
colnames(tmp_vcf_data)[4]="A2"
colnames(tmp_vcf_data)[5]="A1"
colnames(tmp_vcf_data)[10]="mix"
tmp_vcf_data$b=sub("\\:.*", "", tmp_vcf_data$mix)
tmp_vcf_data$SE.1=str_before_nth(tmp_vcf_data$mix, ":", 2)
tmp_vcf_data$se=sub('.*:', '', tmp_vcf_data$SE.1)
tmp_vcf_data$p.1=str_before_nth(tmp_vcf_data$mix, ":", 3)
tmp_vcf_data$p=sub('.*:', '', tmp_vcf_data$p.1)
tmp_vcf_data$eaf.1=str_before_nth(tmp_vcf_data$mix, ":", 4)
tmp_vcf_data$eaf=sub('.*:', '', tmp_vcf_data$eaf.1)
# tmp_vcf_data$N.1=str_before_nth(tmp_vcf_data$mix, ":", 5)
# tmp_vcf_data$N=sub('.*:', '', tmp_vcf_data$N.1)
tmp_vcf_data$p=as.numeric(as.character(tmp_vcf_data$p))
tmp_vcf_data$b=as.numeric(as.character(tmp_vcf_data$b))
#tmp_vcf_data$N=as.numeric(as.character(tmp_vcf_data$N))
tmp_vcf_data$eaf=as.numeric(as.character(tmp_vcf_data$eaf))
tmp_vcf_data$se=as.numeric(as.character(tmp_vcf_data$se))
summary(tmp_vcf_data$p)
tmp_vcf_data$pval=10^(-tmp_vcf_data$p)
tmp_vcf_data$pval2=10^(tmp_vcf_data$p)
summary(tmp_vcf_data$pval)
summary(tmp_vcf_data$p)
tmp_vcf_data$p=tmp_vcf_data$pval
tmp_vcf_data$n=482730
tmp_vcf_data$freq=tmp_vcf_data$eaf
# Select asd variables for SMR
myvars=c("SNP", "A1", "A2", "freq", "b", "se", "p", "n" )
smr_parkinsons=tmp_vcf_data[myvars]
smr_parkinsons=smr_parkinsons[order(smr_parkinsons$p),]
write.table(smr_parkinsons, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/smr_parkinsons.txt",
sep="\t", row.names = FALSE, col.names = TRUE)<!-- sed 's/\"//g' smr_parkinsons.txt > smr_parkinsons -->park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/smr_parkinsons.txt")
park$rownumber = 1:nrow(park)
park=as.data.frame(park)
t.first <- park[match(unique(park$SNP), park$SNP),]
t.first$rownumber=NULL
write.table(t.first, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/park2.txt",
sep="\t", row.names = FALSE, col.names = TRUE)<!-- sed 's/\"//g' park2.txt > smr_park2 -->
# SMR of eqtlGen, GTEx brain, and eMeta brain on Parkinsons
# Run with a bash (slurm / sbatch) script to have the "out" (.o) file assuring completion
# smr_Linux smr --bfile EUR --gwas-summary smr_park2 --beqtl-summary cis-eQTLs-full_eQTLGen_AF_incl_nr_formatted_20191212.new.txt_besd-dense --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_eqtlgen_parkinsons --thread-num 10
# smr_Linux smr --bfile EUR --gwas-summary smr_park2 --beqtl-summary Brain-eMeta --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_Brain_eMeta_parkinsons --thread-num 10
# smr_Linux smr --bfile EUR --gwas-summary smr_park2 --beqtl-summary GTEx_Brain_std --peqtl-smr 5e-4 --out /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_Brain_GTEx_parkinsons --thread-num 10
#cd /n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons
setwd("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons")
park_smr=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_all_smr_parkinsons_SR2_parkOnly2.smr")
# rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
park_smr$tissue=sub('.-*', '-', park_smr$probeID)
park_smr$tissue=str_after_nth(park_smr$probeID, "_", 1)
park_smr=as.data.frame(park_smr)
dim(park_smr)# GPNMB=park_smr[which(park_smr$Gene=="GPNMB"),]# Select the NUC set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_park_smr=park_smr[which(park_smr$Gene %in% combo4$gene),]
rp_park_smr=rp_park_smr[order(rp_park_smr$p_SMR),]
write.table(rp_park_smr, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/rp_park_smr.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)# Read-in emeta brain results
park_emeta=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_Brain_eMeta_parkinsons.smr")
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_park_emeta=park_emeta[which(park_emeta$Gene %in% combo4$gene),]
rp_park_emeta=rp_park_emeta[order(rp_park_emeta$p_SMR),]
rp_park_emeta=as.data.frame(rp_park_emeta)
rp_park_emeta$tissue="eMeta_(cross-brain)"
write.table(rp_park_emeta, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/rp_park_emeta.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
# Read-in GTEx brain results
park_GTEx=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_Brain_GTEx_parkinsons.smr")
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_park_GTEx=park_GTEx[which(park_GTEx$Gene %in% combo4$gene),]
rp_park_GTEx=rp_park_GTEx[order(rp_park_GTEx$p_SMR),]
rp_park_GTEx=as.data.frame(rp_park_GTEx)
rp_park_GTEx$tissue="GTEx_(cross-brain)"
write.table(rp_park_GTEx, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/rp_park_GTEx.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
# Read-in eqtlgen results
park_eqtlgen=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/mysmr_eqtlgen_parkinsons.smr")
# Get the gene symbols
ensemblsIDS=park_eqtlgen$Gene
park_eqtlgen$GENEID=park_eqtlgen$Gene
# Remove the version numbers at the end of the ensembl names
#gtex_pairs_rsid$GENEID <- gsub('\\..+$', '', gtex_pairs_rsid$phenotype_id)
# Fetch the gene names: convert from ensembl.gene to gene.symbol
ensembl.genes=as.character(park_eqtlgen$GENEID)
geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype =
"GENEID", columns = c("SYMBOL","GENEID"))
# Merge in the gene names
park_eqtlgen2 <- merge(x = geneIDs1,
y = park_eqtlgen,
by.x="GENEID",
by.y="GENEID")
park_eqtlgen2$GENEID=NULL
park_eqtlgen2$Gene=park_eqtlgen2$SYMBOL
park_eqtlgen2$SYMBOL=NULL
head(park_eqtlgen)
park_eqtlgen=park_eqtlgen2
# Select the RP set
combo4=read.csv('/n/home04/cdadams/MR-eqtlGen-long/combo4.csv')
rp_park_eqtlgen=park_eqtlgen[which(park_eqtlgen$Gene %in% combo4$gene),]
rp_park_eqtlgen=rp_park_eqtlgen[order(rp_park_eqtlgen$p_SMR),]
rp_park_eqtlgen=as.data.frame(rp_park_eqtlgen)
rp_park_eqtlgen$tissue="Whole_Blood_(eqtlGen)"
write.table(rp_park_eqtlgen, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/rp_park_eqtlgen.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
WB_park=rbind(rp_park_smr,rp_park_emeta,rp_park_GTEx,rp_park_eqtlgen)
head(WB_park)
dim(WB_park)
WB_park$GWAS="Parkinsons"
write.table(WB_park, file ="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/RP_WB_park.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)WB_asd=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/RP_WB_asd.txt")
WB_alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/RP_WB_alz.txt")
WB_park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/RP_WB_park.txt")
three=rbind(WB_asd,WB_alz,WB_park)
table(three$GWAS)
write.table(three, file="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/SMR_three.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
three=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/SMR_three.txt")p=three$p_SMR
bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
three$bonferroni_3outcomes=bonferroni
three_bons=three[which(three$bonferroni_3outcomes<0.05),]
three_p_sigs=three[which(three$p_SMR<0.05),]
three_p_sigs=three_p_sigs[order(three_p_sigs$p_SMR),]
p=three$p_SMR
fdr_3outcomes=p.adjust(p, method = "fdr", n = length(p))
three$fdr_3outcomes=fdr_3outcomes
three_fdr_sigs=three[which(three$fdr_3outcomes<0.05),]
three=as.data.frame(three)three$SNP=three$topSNP
three$beta.outcome=three$b_GWAS
three$se.outcome=three$se_GWAS
three$effect_allele.outcome=three$A1
three$other_allele.outcome=three$A2
three$eaf.outcome=three$Freq
three$outcome="_"
three$id.outcome="_"
three$beta.exposure=three$b_eQTL
three$se.exposure=three$se_eQTL
three$effect_allele.exposure=three$A1
three$other_allele.exposure=three$A2
three$eaf.exposure=three$Freq
three$exposure=three$Gene
three$id.exposure=three$Gene
three$mr_keep="TRUE"
three$mr_keep=as.logical(three$mr_keep)
# Clean up the tissue names for "three" -----------------------------------
# Make the tissue names nice:
table(three$tissue)
# firstup <- function(x) {
# substr(x, 1, 1) <- toupper(substr(x, 1, 1))
# x
# }
# three$tissue=firstup(three$tissue)
three$tissue2=str_remove(three$tissue, "Brain_")
three$tissue2=gsub("_", " ", three$tissue2)
three$tissue2=tolower(three$tissue2)
table(three$tissue2)
three <- transform(three,tissue2=plyr::revalue(tissue2,c("adipose subcutaneous"="adipose (subcutaneous)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("adipose visceral omentum"="adipose (visceral omentum)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("anterior cingulate cortex ba24"="anterior cingulate cortex (ba24)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("artery aorta"="artery (aorta)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("artery tibial"="artery (tibial)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("breast mammary tissue"="breast (mammary tissue)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("caudate basal ganglia"="caudate (basal ganglia)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("putamen basal ganglia"="putamen (basal ganglia)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("cells cultured fibroblasts"="cells (cultured fibroblasts)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("cells ebv-transformed lymphocytes"="cells (ebv-transformed lymphocytes)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("whole blood (gtex"="whole blood (gtex)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("nucleus accumbens basal ganglia"="nucleus accumbens (basal ganglia)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("frontal cortex ba9"="frontal cortex (ba9)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("colon sigmoid"="colon (sigmoid)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("colon transverse"="colon (transverse)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("esophagus gastroesophageal junction"="esophagus (gastroesophageal junction)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("esophagus mucosa"="esophagus (mucosa)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("esophagus muscularis"="esophagus (muscularis)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("heart atrial appendage"="heart (atrial appendage)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("heart left ventricle"="heart (left ventricle)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("skin not sun exposed suprapubic"="skin (not sun exposed suprapubic)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("skin sun exposed lower leg"="skin (sun exposed lower leg)")))
three <- transform(three,tissue2=plyr::revalue(tissue2,c("spinal cord cervical c-1"="spinal cord cervica (c-1)")))
table(three$tissue2)# ASD
ASD=three[three$GWAS=="ASD",]
p=ASD$p_SMR
fdr_ASD=p.adjust(p, method = "fdr", n = length(p))
ASD$fdr_ASD=fdr_ASD
ASD_fdr_set=ASD[which(ASD$fdr_ASD<0.05),]
dim(ASD_fdr_set)
ASD_fdr_set
# Alzheimer's
alz=three[three$GWAS=="Alzheimers",]
p=alz$p_SMR
fdr_alz=p.adjust(p, method = "fdr", n = length(p))
alz$fdr_alz=fdr_alz
alz_fdr_set=alz[which(alz$fdr_ASD<0.05),]
dim(alz_fdr_set)
alz_fdr_set
summary(alz$p_SMR)
# Parkinson's
park=three[three$GWAS=="Parkinsons",]
p=park$p_SMR
fdr_park=p.adjust(p, method = "fdr", n = length(p))
park$fdr_park=fdr_park
park_fdr_set=park[which(park$fdr_park<0.05),]
dim(park_fdr_set)
park_fdr_set
# MALSU1=park[park$Gene=="MALSU1",]
# summary(MALSU1$b_SMR)
# RPS23=three[three$Gene=="RPS23",]
# summary(RPS23$b_SMR)
write.table(three, file="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/SMR_three.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
write.table(park, file="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/park_for_ivw.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
write.table(alz, file="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/alz_for_ivw.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)
write.table(ASD, file="/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/ASD_for_ivw_1030.txt",
sep = "\t",
row.names = FALSE, col.names = TRUE, quote=FALSE)rm(list = ls())
ASD=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/ASD_for_ivw_1030.txt")
ASD$outcome="Autism spectrum disorder (ASD)"
ASD$id.outcome="Autism spectrum disorder (ASD)"
# ASD$tissue3=tolower(ASD$tissue)
# ASD_brain <- ASD[grep("brain", ASD$tissue3), ]
# ASD_brain=as.data.frame(ASD_brain)
ASD_brain <- ASD[grep("Brain", ASD$tissue), ]
ASD_brain=as.data.frame(ASD_brain)
ASD_heidi=ASD_brain[which(ASD_brain$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('ASD_heidi'))])
split <- split(ASD_heidi, ASD_heidi$Gene)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(ASD_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/by_gene/",
names(files[i]), "_ASD_IVW_results_heidi_in_brain_by_gene.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
df_ivw_bons=df_ivw[which(df_ivw$p_bonferroni<0.05),]
table(df_ivw_bons$id.exposure)
# Save the bonferroni IVW results
write.csv(df_ivw_bons, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/by_gene/df_ivw_ASD_heidi_brain_by_gene_bonferroni.csv")
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/by_gene/df_ivw_ASD_heidi_brain_by_gene.csv")
df_ivw_ASD_heidi_brain=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/by_gene/df_ivw_ASD_heidi_brain_by_gene.csv")rm(list = ls())
ASD=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/ASD_for_ivw_1030.txt")
ASD$outcome="Autism spectrum disorder (ASD)"
ASD$id.outcome="Autism spectrum disorder (ASD)"
# ASD$tissue3=tolower(ASD$tissue)
# ASD_brain <- ASD[grep("brain", ASD$tissue3), ]
# ASD_brain=as.data.frame(ASD_brain)
ASD_brain <- ASD[grep("Brain", ASD$tissue), ]
ASD_brain=as.data.frame(ASD_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
ASD_brain_rps_mrps=ASD_brain[which(ASD_brain$Gene%in%rp_list_expanded$all_RP),]
ASD_brain_rps=subset(ASD_brain_rps_mrps,grepl("^RP",Gene,ignore.case = T))
ASD_heidi=ASD_brain_rps[which(ASD_brain_rps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
ASD_heidi$id.exposure="RPs"
ASD_heidi$exposure=ASD_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('ASD_heidi'))])
split <- split(ASD_heidi, ASD_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(ASD_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_rps/",
names(files[i]), "_ASD_IVW_results_by_tissue_RPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_rps/df_ivw_ASD_heidi_brain_RPsonly.csv")rm(list = ls())
ASD=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/ASD_for_ivw_1030.txt")
ASD$outcome="Autism spectrum disorder (ASD)"
ASD$id.outcome="Autism spectrum disorder (ASD)"
# ASD$tissue3=tolower(ASD$tissue)
# ASD_brain <- ASD[grep("brain", ASD$tissue3), ]
# ASD_brain=as.data.frame(ASD_brain)
ASD_brain <- ASD[grep("Brain", ASD$tissue), ]
ASD_brain=as.data.frame(ASD_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
ASD_brain_rps_mrps=ASD_brain[which(ASD_brain$Gene%in%rp_list_expanded$all_RP),]
ASD_brain_mrps=subset(ASD_brain_rps_mrps,grepl("^MRP",Gene,ignore.case = T))
ASD_heidi=ASD_brain_mrps[which(ASD_brain_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
ASD_heidi$id.exposure="MRPs"
ASD_heidi$exposure=ASD_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('ASD_heidi'))])
split <- split(ASD_heidi, ASD_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(ASD_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_mrps/",
names(files[i]), "_ASD_IVW_results_by_tissue_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_mrps/df_ivw_ASD_heidi_brain_MRPsonly.csv")rm(list = ls())
ASD=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/ASD_for_ivw_1030.txt")
ASD$outcome="Autism spectrum disorder (ASD)"
ASD$id.outcome="Autism spectrum disorder (ASD)"
# ASD$tissue3=tolower(ASD$tissue)
# ASD_brain <- ASD[grep("brain", ASD$tissue3), ]
# ASD_brain=as.data.frame(ASD_brain)
ASD_brain <- ASD[grep("Brain", ASD$tissue), ]
ASD_brain=as.data.frame(ASD_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
ASD_brain_rps_mrps=ASD_brain[which(ASD_brain$Gene%in%rp_list_expanded$all_RP),]
ASD_heidi=ASD_brain_rps_mrps[which(ASD_brain_rps_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
ASD_heidi$id.exposure="RPs and MRPs"
ASD_heidi$exposure=ASD_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('ASD_heidi'))])
split <- split(ASD_heidi, ASD_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(ASD_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_rps_mrps/",
names(files[i]), "_ASD_IVW_results_by_tissue_RPs_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR correction
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ASD_IVW/across_rps_mrps/df_ivw_ASD_heidi_brain_RPs_MRPs.csv")rm(list = ls())
alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/alz_for_ivw.txt")
alz$outcome="Alzheimer's disease"
alz$id.outcome="Alzheimer's disease"
# alz$tissue3=tolower(alz$tissue)
# alz_brain <- alz[grep("brain", alz$tissue3), ]
# alz_brain=as.data.frame(alz_brain)
alz_brain <- alz[grep("Brain", alz$tissue), ]
alz_brain=as.data.frame(alz_brain)
alz_heidi=alz_brain[which(alz_brain$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('alz_heidi'))])
split <- split(alz_heidi, alz_heidi$Gene)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(alz_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/by_gene/",
names(files[i]), "_alz_IVW_alz_results_heidi_in_brain_by_gene.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
df_ivw_bons=df_ivw[which(df_ivw$p_bonferroni<0.05),]
table(df_ivw_bons$id.exposure)
# Save the bonferroni IVW results
write.csv(df_ivw_bons, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/by_gene/df_ivw_alz_heidi_brain_by_gene_bonferroni.csv")
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/by_gene/df_ivw_ASD_heidi_brain_by_gene.csv")
df_ivw_ASD_heidi_brain=read.csv("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/by_gene/df_ivw_ASD_heidi_brain_by_gene.csv")rm(list = ls())
alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/alz_for_ivw.txt")
alz$outcome="Alzheimer's disease"
alz$id.outcome="Alzheimer's disease"
# alz$tissue3=tolower(alz$tissue)
# alz_brain <- alz[grep("brain", alz$tissue3), ]
# alz_brain=as.data.frame(alz_brain)
alz_brain <- alz[grep("Brain", alz$tissue), ]
alz_brain=as.data.frame(alz_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
alz_brain_rps_mrps=alz_brain[which(alz_brain$Gene%in%rp_list_expanded$all_RP),]
alz_brain_rps=subset(alz_brain_rps_mrps,grepl("^RP",Gene,ignore.case = T))
alz_heidi=alz_brain_rps[which(alz_brain_rps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
alz_heidi$id.exposure="RPs"
alz_heidi$exposure=alz_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('alz_heidi'))])
split <- split(alz_heidi, alz_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(alz_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_rps/",
names(files[i]), "_alz_IVW_results_by_tissue_RPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_rps/df_ivw_alz_heidi_brain_RPsonly.csv")rm(list = ls())
alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/alz_for_ivw.txt")
alz$outcome="Alzheimer's disease"
alz$id.outcome="Alzheimer's disease"
# alz$tissue3=tolower(alz$tissue)
# alz_brain <- alz[grep("brain", alz$tissue3), ]
# alz_brain=as.data.frame(alz_brain)
alz_brain <- alz[grep("Brain", alz$tissue), ]
alz_brain=as.data.frame(alz_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
alz_brain_rps_mrps=alz_brain[which(alz_brain$Gene%in%rp_list_expanded$all_RP),]
alz_brain_mrps=subset(alz_brain_rps_mrps,grepl("^MRP",Gene,ignore.case = T))
alz_heidi=alz_brain_mrps[which(alz_brain_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
alz_heidi$id.exposure="MRPs"
alz_heidi$exposure=alz_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('alz_heidi'))])
split <- split(alz_heidi, alz_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(alz_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_mrps/",
names(files[i]), "_alz_IVW_results_by_tissue_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_mrps/df_ivw_ALZ_heidi_brain_MRPsonly.csv")rm(list = ls())
alz=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/alzheimers/alz_for_ivw.txt")
alz$outcome="Alzheimer's disease"
alz$id.outcome="Alzheimer's disease"
# alz$tissue3=tolower(alz$tissue)
# alz_brain <- alz[grep("brain", alz$tissue3), ]
# alz_brain=as.data.frame(alz_brain)
alz_brain <- alz[grep("Brain", alz$tissue), ]
alz_brain=as.data.frame(alz_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
alz_brain_rps_mrps=alz_brain[which(alz_brain$Gene%in%rp_list_expanded$all_RP),]
alz_heidi=alz_brain_rps_mrps[which(alz_brain_rps_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
alz_heidi$id.exposure="RPs and MRPs"
alz_heidi$exposure=alz_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('alz_heidi'))])
split <- split(alz_heidi, alz_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(alz_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_rps_mrps/",
names(files[i]), "_alz_IVW_results_by_tissue_RPs_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR correction
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_ALZ_IVW/across_rps_mrps/df_ivw_alz_heidi_brain_RPs_MRPs.csv")rm(list = ls())
park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/park_for_ivw.txt")
park$outcome="Parkinson's disease"
park$id.outcome="Parkinson's disease"
# park$tissue3=tolower(park$tissue)
# park_brain <- park[grep("brain", park$tissue3), ]
# park_brain=as.data.frame(park_brain)
park_brain <- park[grep("Brain", park$tissue), ]
park_brain=as.data.frame(park_brain)
park_heidi=park_brain[which(park_brain$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('park_heidi'))])
split <- split(park_heidi, park_heidi$Gene)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(park_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/by_gene/",
names(files[i]), "_park_IVW_park_results_heidi_in_brain_by_gene.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
df_ivw_bons=df_ivw[which(df_ivw$p_bonferroni<0.05),]
table(df_ivw_bons$id.exposure)
# Save the bonferroni IVW results
write.csv(df_ivw_bons, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/by_gene/df_ivw_park_heidi_brain_by_gene_bonferroni.csv")
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/by_gene/df_ivw_park_heidi_brain_by_gene.csv")
df_ivw_park_heidi_brain=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/by_gene/df_ivw_park_heidi_brain_by_gene.csv")rm(list = ls())
park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/park_for_ivw.txt")
park$outcome="Parkinson's disease"
park$id.outcome="Parkinson's disease"
# park$tissue3=tolower(park$tissue)
# park_brain <- park[grep("brain", park$tissue3), ]
# park_brain=as.data.frame(park_brain)
park_brain <- park[grep("Brain", park$tissue), ]
park_brain=as.data.frame(park_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
park_brain_rps_mrps=park_brain[which(park_brain$Gene%in%rp_list_expanded$all_RP),]
park_brain_rps=subset(park_brain_rps_mrps,grepl("^RP",Gene,ignore.case = T))
park_heidi=park_brain_rps[which(park_brain_rps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
park_heidi$id.exposure="RPs"
park_heidi$exposure=park_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('park_heidi'))])
split <- split(park_heidi, park_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(park_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps/",
names(files[i]), "_Park_IVW_results_by_tissue_RPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps/df_ivw_park_heidi_brain_RPsonly.csv")
df_ivw_park_heidi_brain_RPsonly=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps/df_ivw_park_heidi_brain_RPsonly.csv")
df_ivw_park_heidi_brain_RPsonlyrm(list = ls())
park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/park_for_ivw.txt")
park$outcome="Parkinson's disease"
park$id.outcome="Parkinson's disease"
# park$tissue3=tolower(park$tissue)
# park_brain <- park[grep("brain", park$tissue3), ]
# park_brain=as.data.frame(park_brain)
park_brain <- park[grep("Brain", park$tissue), ]
park_brain=as.data.frame(park_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
park_brain_rps_mrps=park_brain[which(park_brain$Gene%in%rp_list_expanded$all_RP),]
park_brain_mrps=subset(park_brain_rps_mrps,grepl("^MRP",Gene,ignore.case = T))
park_heidi=park_brain_mrps[which(park_brain_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
park_heidi$id.exposure="MRPs"
park_heidi$exposure=park_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('park_heidi'))])
split <- split(park_heidi, park_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(park_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_mrps/",
names(files[i]), "_park_IVW_results_by_tissue_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_mrps/df_ivw_park_heidi_brain_MRPsonly.csv")
df_ivw_park_heidi_brain_MRPsonly=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_mrps/df_ivw_park_heidi_brain_MRPsonly.csv")
df_ivw_park_heidi_brain_MRPsonlyrm(list = ls())
park=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/parkinsons/park_for_ivw.txt")
park$outcome="Parkinson's disease"
park$id.outcome="Parkinson's disease"
# park$tissue3=tolower(park$tissue)
# park_brain <- park[grep("brain", park$tissue3), ]
# park_brain=as.data.frame(park_brain)
park_brain <- park[grep("Brain", park$tissue), ]
park_brain=as.data.frame(park_brain)
rp_list_expanded=read.csv('/n/home04/cdadams/ld/ldsc/expanded_verified.csv')
park_brain_rps_mrps=park_brain[which(park_brain$Gene%in%rp_list_expanded$all_RP),]
park_heidi=park_brain_rps_mrps[which(park_brain_rps_mrps$p_HEIDI>0.05),]#& three_for_IVW$p_SMR<0.05
park_heidi$id.exposure="RPs and MRPs"
park_heidi$exposure=park_heidi$tissue2
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('park_heidi'))])
split <- split(park_heidi, park_heidi$tissue2)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(park_heidi)
files <- mget(ls())
length(files)
for (i in 1:length(files)){
write.csv(files[[i]],
paste("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps_mrps/",
names(files[i]), "_park_IVW_results_by_tissue_RPs_MRPs.csv", sep = ""))
}
result <- list()
for(i in 1:length(files))
{
result[[i]] <- mr(files[[i]])
}
result[1]
df <- do.call("rbind", result)
df_ivw=df[which(df$method=="Inverse variance weighted"),]
# Bonferroni correction
p=df_ivw$pval
p_bonferroni=p.adjust(p, method = "bonferroni", n = length(p))
df_ivw$p_bonferroni=p_bonferroni
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# FDR correction
p=df_ivw$pval
p_fdr=p.adjust(p, method = "fdr", n = length(p))
df_ivw$p_fdr=p_fdr
df_ivw=df_ivw[order(df_ivw$p_bonferroni),]
# Save the full IVW results
write.csv(df_ivw, "/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps_mrps/df_ivw_park_heidi_brain_RPs_MRPs.csv")
df_ivw_park_heidi_brain_RPs_MRPs=fread("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/Brain_Park_IVW/across_rps_mrps/df_ivw_park_heidi_brain_RPs_MRPs.csv")
df_ivw_park_heidi_brain_RPs_MRPsivws=read.csv("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/full_ivw_res.csv")
rps23_ivw=ivws[which(ivws$id.exposure=="RPS23"),]
df_add_ivw=read.csv("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/WB_for_IVW_1026.csv")
RPS23_ivw=df_ivw[which(df_ivw$id.exposure=="RPS23"),]
RPS23=RPS23_ivw %>% add_row(tissue = "Inverse variance weighted (IVW) meta-analysis", topSNP="IVW",
b_SMR = 0.05982229, se_SMR=0.003852278, p_SMR=6.947896e-52)
forest=forestplot(
df = RPS23,
name = tissue,
estimate = b_SMR,
se = se_SMR,
pvalue = p_SMR,
#psignif = 3.608597e-07,
colour = topSNP, #Trait
#shape = id.exposure,
xlab = "SMR results: betas and 95% confidence intervals per tissue and \n IVW meta-analysis (IVW Bonferroni-corrected P-value=6.947896e-52)",
title = RPS23$Gene,
logodds = FALSE,
# xlim = c(0.6, 1.4),
# You can explicitly define x-tick breaks
# xtickbreaks = c(0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4)
)
forest
ggsave("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/full_body_forest_RPS23.png",
units="in",width=8, height=8, dpi=300)setwd("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res")
WB_for_IVW=read.csv("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/WB_for_IVW_1026.csv")
WB_heidi=WB_for_IVW[which(WB_for_IVW$p_HEIDI>0.05 ),]#& WB_for_IVW$p_SMR<0.05
# Remove all objects but ...
rm(list= ls()[!(ls() %in% c('WB_heidi'))])
split <- split(WB_heidi, WB_heidi$Gene)
list2env(split, envir= .GlobalEnv) #split the list into separate datasets
rm(split)
rm(WB_heidi)
files <- mget(ls())
length(files)
names(files)
# Remove the empty data frames
# Create a function that returns a logical value
isEmpty <- function(x) {
is.data.frame(x) && nrow(x) == 0L
}
# Apply it over the environment
empty <- unlist(eapply(.GlobalEnv, isEmpty))
# Remove the empties
rm(list = names(empty)[empty])
rm(empty)
rm(isEmpty)
rm(files)
# use "zz" as a work around for forestplot recognizing it isn't a df
zz <- mget(ls())
length(zz)
plot_list = list()
for(i in 1:length(zz)){
p=forestplot(
df = zz[[i]],
name = tissue2,
estimate = b_SMR,
se = se_SMR,
pvalue = p_SMR,
#psignif = 3.608597e-07,
colour = SNP,
xlab = "MR results: betas and 95% confidence intervals",
title = zz[[i]]$Gene,
logodds = FALSE
)
plot_list[[i]] = p
print(plot_list[[i]])
ggsave(plot_list[[i]],
file=paste0("/n/holylfs/LABS/lemos_lab/Users/cdadams/gtex_EUR/all_tissue_SMR/job_res/forest_all_body/",
names(zz)[[i]],".png"),
units="in",width=8, height=15, dpi=300)
}
# combo_density <- ggarrange(life, health, nrow = 1, ncol = 2, legend = NULL, legend.grob = NULL)
# combo_density