In [2]:
library(biomaRt) # Load the biomaRt library for accessing biological data through BioMart
library(data.table) # Load the data.table library for efficient data manipulation
library(ggplot2) # Load the ggplot2 library for creating data visualizations
library(readxl) # Load the readxl library for reading Excel files into R
In [3]:
setwd("/scratch/silo1/rameez/ldscore/Geneweaver")
GeneWeaver database¶
In [4]:
# Read an Excel file named "Withdrawal Human Orthologs.xlsx" from the specified path
Jason_list <- read_excel("/projects/bga_lab/DATA_REPOSITORIES/COGA_WITHDRAWAL/UM_IMP_RESULTS/MERGED_CLEAN/SmokeScreen/Withdrawal Human Orthologs.xlsx",sheet = 1, skip = 7)
# Convert the read data into a data frame
Jason_list <- as.data.frame(Jason_list)
# Display the first few rows and dimensions of the data frame
head(Jason_list); dim(Jason_list)
# Create a new column named "GENE" containing lowercase versions of values from "Column 2"
Jason_list$GENE <- tolower(Jason_list$'Column 2')
# Extract columns "GENE" and "Hippocampus Alcohol Withdrawal" into a new data frame
J_list_withdrawal <- Jason_list[,c("GENE","Hippocampus Alcohol Withdrawal")]
# Filter rows where the "Hippocampus Alcohol Withdrawal" column is equal to 1
J_list_withdrawal <- J_list_withdrawal[J_list_withdrawal$'Hippocampus Alcohol Withdrawal' == 1,]
# Display the first few rows and dimensions of the filtered data frame
J_list_withdrawal$GENE <- toupper(J_list_withdrawal$GENE )
head(J_list_withdrawal); dim(J_list_withdrawal)
Option #1: Utilize the BioMart packages in R to extract chromosome, start/end positions.¶
In [6]:
# Define the biomart connection for GRCh37
grch37 = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org", path="/biomart/martservice", dataset="hsapiens_gene_ensembl")
# Retrieve gene information from GRCh37 biomart
gene_37 <-getBM(attributes=c('ensembl_gene_id',"ensembl_gene_id_version","hgnc_symbol",'chromosome_name','start_position','end_position'), mart = grch37)
# Filter genes located on chromosomes 1 to 22
gene_37 <- gene_37[gene_37$chromosome_name == 1 | gene_37$chromosome_name == 2 | gene_37$chromosome_name == 3 | gene_37$chromosome_name == 4 |
gene_37$chromosome_name == 5 | gene_37$chromosome_name == 6 | gene_37$chromosome_name == 7 | gene_37$chromosome_name == 8 |
gene_37$chromosome_name == 9 | gene_37$chromosome_name == 10 | gene_37$chromosome_name == 11 | gene_37$chromosome_name == 12 |
gene_37$chromosome_name == 13 | gene_37$chromosome_name == 14 | gene_37$chromosome_name == 15 | gene_37$chromosome_name == 16 |
gene_37$chromosome_name == 17 | gene_37$chromosome_name == 18 | gene_37$chromosome_name == 19 | gene_37$chromosome_name == 20 |
gene_37$chromosome_name == 21 | gene_37$chromosome_name == 22 ,]
# Convert the result to a data frame
gene_37 <- as.data.frame(gene_37)
# Display the first few rows and dimensions of the data frame
head(gene_37); dim(gene_37)
In [11]:
# Merge the J_list_withdrawal and gene_37 data frames by the 'GENE' column without sorting
gene_37_J_list_withdrawal <- merge(J_list_withdrawal,gene_37,by.x = "GENE", by.y = "hgnc_symbol", sort = FALSE)
# Display the first few rows and dimensions of the merged data frame
head(gene_37_J_list_withdrawal); dim(gene_37_J_list_withdrawal)
In [14]:
# Select specific columns 'chromosome_name', 'start_position', and 'end_position' from gene_37_J_list_withdrawal
gene_37_J_list_withdrawal_sub <- gene_37_J_list_withdrawal[,c("chromosome_name","start_position","end_position")]
# Prepend "chr" to the 'chromosome_name' column values
gene_37_J_list_withdrawal_sub$chromosome_name <- paste("chr",gene_37_J_list_withdrawal_sub$chromosome_name,sep = "")
# Rename the column names to 'CHR', 'START', and 'END'
colnames(gene_37_J_list_withdrawal_sub) <- c("CHR","START","END")
# Display the first few rows and dimensions of the subsetted data frame
head(gene_37_J_list_withdrawal_sub); dim(gene_37_J_list_withdrawal_sub)
# Write the subsetted data frame to a tab-separated text file named "J_list_withdrawal_geneset.txt"
fwrite(gene_37_J_list_withdrawal_sub,"J_list_withdrawal_geneset.txt",sep = "\t")
Option #2: Use the following link (https://biit.cs.ut.ee/gprofiler/convert) to extract the ENSC format.¶
In [15]:
fwrite(as.data.frame(J_list_withdrawal[,1]),"J_list_withdrawal_gene.txt", col.names = F)
In [16]:
IRdisplay::display_png(file = "gene_list.png")
In [24]:
J_list_withdrawal_ENSG <- fread("gProfiler_hsapiens_7-24-2023_12-41-16 PM.csv")
head(J_list_withdrawal_ENSG); dim(J_list_withdrawal_ENSG)
sum(is.na(J_list_withdrawal_ENSG$converted_alias))
sum(J_list_withdrawal_ENSG$converted_alias == "None")
fwrite(as.data.frame(J_list_withdrawal_ENSG[,2]),"J_list_withdrawal_ENSG.txt", col.names = F)
Annot¶
In [35]:
# Run In Linux
mkdir J_list_withdrawal
for i in {1..22}
do
/projects/bga_lab/Rameez/Tool/miniconda3/envs/ldsc/bin/python \
/projects/bga_lab/Rameez/Tool/ldsc/make_annot.py \
--gene-set-file /scratch/silo1/rameez/ldscore/Geneweaver/J_list_withdrawal_ENSG.txt \
--gene-coord-file /scratch/silo1/rameez/ldscore/analysis/reference/make_annot_sample_files/ENSG_coord.txt \
--windowsize 100000 \
--bimfile /scratch/silo1/rameez/ldscore/analysis/reference/1000G_EUR_Phase3_plink/1000G.EUR.QC.${i}.bim \
--annot-file J_list_withdrawal/J_list_withdrawal.${i}.annot.gz
done
Summary Statistics¶
In [38]:
sum_stat_ea <- fread("/projects/bga_lab/DATA_REPOSITORIES/COGA_WITHDRAWAL/UM_IMP_RESULTS/MERGED_CLEAN_V2/update_result/ea_sub/MLMA/COGA_factor_res_age_sex_PC14_cov_cohort_mlma.mlma")
sum_stat_ea <- as.data.frame(sum_stat_ea)
rs <- fread("/projects/bga_lab/DATA_REPOSITORIES/COGA_WITHDRAWAL/UM_IMP_RESULTS/MERGED_CLEAN_V2/update_result/ea_sub/H_magma/sumStat_rs.txt")
rs <- as.data.frame(rs)
sum_stat_ea <- sum_stat_ea[order(sum_stat_ea$p),]
# 7320 Total sample
In [39]:
sum_stat_ea_rs <- merge(sum_stat_ea,rs, by.x = "p", by.y = "P", sort = FALSE)
head(sum_stat_ea_rs); dim(sum_stat_ea_rs)
In [44]:
sum_stat_ea_rs_sub <- sum_stat_ea_rs[,c("SNP.y","A1","A2","Freq","b","se","p")]
sum_stat_ea_rs_sub <- sum_stat_ea_rs_sub[!duplicated(sum_stat_ea_rs_sub$SNP.y),]
colnames(sum_stat_ea_rs_sub)[1:7] <- c("SNP_RS","Allele1","Allele2","Freq.A1","BETA","SE","P")
head(sum_stat_ea_rs_sub); dim(sum_stat_ea_rs_sub); sum(duplicated(sum_stat_ea_rs_sub$SNP_RS))
fwrite(sum_stat_ea_rs_sub,"sum_stat_ea_rs_sub.txt", sep = "\t")
In [ ]:
# Run in Linux
# format summary stats
/projects/bga_lab/Rameez/Tool/miniconda3/envs/ldsc/bin/python \
/projects/bga_lab/Rameez/Tool/ldsc/munge_sumstats.py \
--sumstats /scratch/silo1/rameez/ldscore/Geneweaver/sum_stat_ea_rs_sub.txt \
--merge-alleles /scratch/silo1/rameez/ldscore/analysis/reference/w_hm3.snplist \
--N 7320 --snp SNP_RS --p P --frq Freq.A1 --a1 Allele1 --a2 Allele2 \
--chunksize 50000 \
--out coga_WD_ea_V2
LD score¶
In [ ]:
# Run in linux
# awk '{if ($1!="SNP") {print $1} }' w_hm3.snplist > listHM3.txt
for i in {1..22}
do
/projects/bga_lab/Rameez/Tool/miniconda3/envs/ldsc/bin/python \
/projects/bga_lab/Rameez/Tool/ldsc/ldsc.py \
--l2 \
--bfile /scratch/silo1/rameez/ldscore/analysis/reference/1000G_EUR_Phase3_plink/1000G.EUR.QC.${i} \
--ld-wind-cm 1 \
--annot /scratch/silo1/rameez/ldscore/Geneweaver/J_list_withdrawal/J_list_withdrawal.${i}.annot.gz \
--thin-annot \
--out J_list_withdrawal/J_list_withdrawal.${i} \
--print-snps /scratch/silo1/rameez/ldscore/analysis/reference/listHM3.txt
done
In [ ]:
Enrishment analysis¶
In [ ]:
/projects/bga_lab/Rameez/Tool/miniconda3/envs/ldsc/bin/python \
/projects/bga_lab/Rameez/Tool/ldsc/ldsc.py \
--h2 /scratch/silo1/rameez/ldscore/Geneweaver/coga_WD_ea_V2.sumstats.gz \
--ref-ld-chr /scratch/silo1/rameez/ldscore/Geneweaver/J_list_withdrawal/J_list_withdrawal.,/scratch/silo1/rameez/ldscore/analysis/reference/1000G_EUR_Phase3_baseline/baseline. \
--w-ld-chr /scratch/silo1/rameez/ldscore/analysis/reference/weights_hm3_no_hla/weights. \
--overlap-annot \
--print-coefficients \
--frqfile-chr /scratch/silo1/rameez/ldscore/analysis/reference/1000G_Phase3_frq/1000G.EUR.QC. \
--out J_list_withdrawal/withdrawal
In [51]:
withdrawal <- fread("/scratch/silo1/rameez/ldscore/Geneweaver/J_list_withdrawal/withdrawal.results")
(withdrawal); dim(withdrawal)
In [52]:
brain <- fread("/scratch/silo1/rameez/ldscore/analysis/Brain_dir/BMI_Brain_V3.results")
(brain); dim(brain)
In [ ]:
List of Variable and memory space¶
In [45]:
var_mem <- NULL
var_name <- NULL
for (i in 1:length(ls())){
var_mem[i] <- object.size(get(ls()[i]))
var_name[i] <- ls()[i]
}
In [46]:
var_mem_dat <- cbind(var_name,var_mem)
var_mem_dat <- as.data.frame(var_mem_dat)
var_mem_dat$var_mem <- as.numeric(var_mem_dat$var_mem)
var_mem_dat$megabytes <- (var_mem_dat$var_mem/1000000)
var_mem_dat$gigabytes <- (var_mem_dat$var_mem/1e+9)
(var_mem_dat)
In [50]:
rm(list=ls())
In [62]:
73.302^2 - ( (1-(73.302^2))/(500-2) )
In [64]:
!plink --help
In [ ]:
C:\Users\ramee\Downloads
In [70]:
dat <- fread("/projects/bga_lab/Rameez/Code/results2 (4).results")
head(dat); dim(dat)
In [ ]: