.Rhistory

getwd()
setwd("../")
getwd()
library(pkgdown)
pkgdown::build_site()
pkgdown::build_site()
output.file<-system.file("test", package = "ProxSeek")
output.file
getwd()
demo.P.PIR<-system.file("data", package = "ProxSeek")
demo.P.PIR
output.file<-paste0(output.file,"/output.csv")
result <- ProxSeek_runProxeQTL(SNP_file_name = SNP_file_name,expression_file_name = expression_file_name,snps_location_file_name = snps_location_file_name,gene_location_file_name = gene_location_file_name,covariates_file_name = covariates_file_name, cell.type = "Mon", chr = 19,
output.file = output.file, n.cluster = 10, P_PIR_folder = demo.P.PIR)
output.file<-system.file("test", package = "ProxSeek")
output.file
getwd()
devtools::build()
install.packages("/media/xingewang/JRlab10TB/Project_Space/Proxi_HIC_Validation/Prepare_CellTypePPIR/ProxSeek_1.0.0.tar.gz", repos = NULL, type = "source")
output.file<-system.file("test", package = "ProxSeek")
output.file<-paste0(output.file,"/output.csv")
output.file
output.file<-system.file("data", package = "ProxSeek")
output.file<-paste0(output.file,"/output.csv")
output.file
devtools::build()
pkgdown::build_site()
cell.type<-"Mon"
SNP_file_name = system.file("data", "chr19_geno2.txt", package = "ProxSeek")
expression_file_name = system.file("data", "GE2.txt", package = "ProxSeek")
snps_location_file_name = system.file("data", "chr19_geno_loc2.txt", package = "ProxSeek")
gene_location_file_name = system.file("data", "geneloc_TSS2.txt", package = "ProxSeek")
cell.type<-"Mon"
SNP_file_name = system.file("data", "chr19_geno2.csv", package = "ProxSeek")
expression_file_name = system.file("data", "GE2.csv", package = "ProxSeek")
snps_location_file_name = system.file("data", "chr19_geno_loc2.csv", package = "ProxSeek")
gene_location_file_name = system.file("data", "geneloc_TSS2.csv", package = "ProxSeek")
covariates_file_name = system.file("data", "Covariates.csv", package = "ProxSeek")
chr<-19
chr<-19
##### 0. check chromosome id
chr.i<-chr
chr.arr<-c(seq(1,22), "X", "Y")
if(!chr.i %in% chr.arr){
stop("Chromosome must be between 1-22 or X and Y!!!")
}
##### 1. read p-pir file and get the furthest distance
# cell.type<-"monocyte"
cell.type.file<-paste0("P_PIR_", cell.type,".csv")
cell.type.file<-system.file("data", cell.type.file, package = "ProxSeek")
# Check if the file exists
if (!file.exists(cell.type.file)) {
stop("Cell type P-PIR file not found: ", cell.type.file)
}
# Attempt to read the file
tryCatch({
cell.hic <- fread(cell.type.file)
}, error = function(e) {
stop("Failed to read cell type P-PIR file: ", e$message)
})
# Validate the file content
requiredColumns <- c("P_chr","P_start", "P_end", "P_gene","PIR_chr","PIR_start", "PIR_end", "PIR_gene")
gc()
library(data.table)
library(stringr)
cell.type<-"Mon"
SNP_file_name = system.file("data", "chr19_geno2.csv", package = "ProxSeek")
expression_file_name = system.file("data", "GE2.csv", package = "ProxSeek")
snps_location_file_name = system.file("data", "chr19_geno_loc2.csv", package = "ProxSeek")
gene_location_file_name = system.file("data", "geneloc_TSS2.csv", package = "ProxSeek")
covariates_file_name = system.file("data", "Covariates.csv", package = "ProxSeek")
chr<-19
##### 0. check chromosome id
chr.i<-chr
chr.arr<-c(seq(1,22), "X", "Y")
if(!chr.i %in% chr.arr){
stop("Chromosome must be between 1-22 or X and Y!!!")
}
##### 1. read p-pir file and get the furthest distance
# cell.type<-"monocyte"
cell.type.file<-paste0("P_PIR_", cell.type,".csv")
cell.type.file<-system.file("data", cell.type.file, package = "ProxSeek")
# Check if the file exists
if (!file.exists(cell.type.file)) {
stop("Cell type P-PIR file not found: ", cell.type.file)
}
# Attempt to read the file
tryCatch({
cell.hic <- fread(cell.type.file)
}, error = function(e) {
stop("Failed to read cell type P-PIR file: ", e$message)
})
# Validate the file content
requiredColumns <- c("P_chr","P_start", "P_end", "P_gene","PIR_chr","PIR_start", "PIR_end", "PIR_gene")
if (!all(requiredColumns %in% colnames(cell.hic))) {
stop("Cell type P-PIR file does not contain required columns.")
}
# get the furthest distance
cell.hic.i<-cell.hic %>% filter(P_chr == chr.i & PIR_chr == chr.i)
cell.hic
class(cell.hic)
cell.hic<-as.data.frame(cell.hic)
# Validate the file content
requiredColumns <- c("P_chr","P_start", "P_end", "P_gene","PIR_chr","PIR_start", "PIR_end", "PIR_gene")
if (!all(requiredColumns %in% colnames(cell.hic))) {
stop("Cell type P-PIR file does not contain required columns.")
}
# get the furthest distance
cell.hic.i<-cell.hic %>% filter(P_chr == chr.i & PIR_chr == chr.i)
dist.df<-ddply(cell.hic.i, .(P_gene), function(df){
### get max distance to TSS of the gene
### down stream
down.idx<-which(df$PIR_start - df$P_end >0)
up.idx<-which(df$PIR_end - df$P_start<0)
if(length(down.idx)!=0){
down.dist<-df$PIR_end - df$P_start
down.dist<-down.dist[down.idx]
max.down.dist<-max(down.dist)
}else{max.down.dist<-0}
if(length(up.idx)!=0){
up.dist<-df$P_end - df$PIR_start
up.dist<-up.dist[up.idx]
max.up.dist<-max(up.dist)
}else{max.up.dist<-0}
return(max(max.down.dist,max.up.dist))
})
library(dplyr)
library(plyr)
cell.type<-"Mon"
SNP_file_name = system.file("data", "chr19_geno2.csv", package = "ProxSeek")
expression_file_name = system.file("data", "GE2.csv", package = "ProxSeek")
snps_location_file_name = system.file("data", "chr19_geno_loc2.csv", package = "ProxSeek")
gene_location_file_name = system.file("data", "geneloc_TSS2.csv", package = "ProxSeek")
covariates_file_name = system.file("data", "Covariates.csv", package = "ProxSeek")
chr<-19
##### 0. check chromosome id
chr.i<-chr
chr.arr<-c(seq(1,22), "X", "Y")
if(!chr.i %in% chr.arr){
stop("Chromosome must be between 1-22 or X and Y!!!")
}
##### 1. read p-pir file and get the furthest distance
# cell.type<-"monocyte"
cell.type.file<-paste0("P_PIR_", cell.type,".csv")
cell.type.file<-system.file("data", cell.type.file, package = "ProxSeek")
# Check if the file exists
if (!file.exists(cell.type.file)) {
stop("Cell type P-PIR file not found: ", cell.type.file)
}
# Attempt to read the file
tryCatch({
cell.hic <- fread(cell.type.file)
}, error = function(e) {
stop("Failed to read cell type P-PIR file: ", e$message)
})
#cell.hic<-as.data.frame(cell.hic)
# Validate the file content
requiredColumns <- c("P_chr","P_start", "P_end", "P_gene","PIR_chr","PIR_start", "PIR_end", "PIR_gene")
if (!all(requiredColumns %in% colnames(cell.hic))) {
stop("Cell type P-PIR file does not contain required columns.")
}
# get the furthest distance
cell.hic.i<-cell.hic %>% filter(P_chr == chr.i & PIR_chr == chr.i)
dist.df<-ddply(cell.hic.i, .(P_gene), function(df){
### get max distance to TSS of the gene
### down stream
down.idx<-which(df$PIR_start - df$P_end >0)
up.idx<-which(df$PIR_end - df$P_start<0)
if(length(down.idx)!=0){
down.dist<-df$PIR_end - df$P_start
down.dist<-down.dist[down.idx]
max.down.dist<-max(down.dist)
}else{max.down.dist<-0}
if(length(up.idx)!=0){
up.dist<-df$P_end - df$PIR_start
up.dist<-up.dist[up.idx]
max.up.dist<-max(up.dist)
}else{max.up.dist<-0}
return(max(max.down.dist,max.up.dist))
})
max.cis.dist<-max(dist.df$V1)
max.cis.dist
library(MatrixEQTL)
##### 2. run MatrixeQTL
# Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
# Genotype file name
SNP_file_name = SNP_file_name
snps_location_file_name = snps_location_file_name
# Gene expression file name
expression_file_name = expression_file_name
gene_location_file_name = gene_location_file_name
# Covariates file name
# Set to character() for no covariates
covariates_file_name = covariates_file_name
# Output file name
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();
# Only associations significant at this level will be saved
pvOutputThreshold_cis = 0.05;
pvOutputThreshold_tra = 0;
# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric();
# Distance for local gene-SNP pairs
cisDist = max.cis.dist;
snps = SlicedData$new();
snps$fileDelimiter = "\t";      # the TAB character
snps$fileOmitCharacters = "NA"; # denote missing values;
snps$fileSkipRows = 1;          # one row of column labels
snps$fileSkipColumns = 1;       # one column of row labels
snps$fileSliceSize = 2000;      # read file in slices of 2,000 rows
snps$LoadFile(SNP_file_name);
gene = SlicedData$new();
gene$fileDelimiter = "\t";      # the TAB character
gene$fileOmitCharacters = "NA"; # denote missing values;
gene$fileSkipRows = 1;          # one row of column labels
gene$fileSkipColumns = 1;       # one column of row labels
gene$fileSliceSize = 2000;      # read file in slices of 2,000 rows
gene$LoadFile(expression_file_name);
cvrt = SlicedData$new();
cvrt$fileDelimiter = "\t";      # the TAB character
cvrt$fileOmitCharacters = "NA"; # denote missing values;
cvrt$fileSkipRows = 1;          # one row of column labels
cvrt$fileSkipColumns = 1;       # one column of row labels
if(length(covariates_file_name)>0) {
cvrt$LoadFile(covariates_file_name);
}
## Run the analysis
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);
me = Matrix_eQTL_main(
snps = snps,
gene = gene,
cvrt = cvrt,
output_file_name     = output_file_name_tra,
pvOutputThreshold     = pvOutputThreshold_tra,
useModel = useModel,
errorCovariance = errorCovariance,
verbose = TRUE,
output_file_name.cis = output_file_name_cis,
pvOutputThreshold.cis = pvOutputThreshold_cis,
snpspos = snpspos,
genepos = genepos,
cisDist = cisDist,
pvalue.hist = "qqplot",
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
unlink(output_file_name_tra);
unlink(output_file_name_cis);
### add required column names c("snps", "gene","Symbol","chr","pos")
snp.loc.df<-fread(snps_location_file_name)
result.df<-me$cis$eqtls
result.df$chr<-snp.loc.df$chr[match(result.df$snps, snp.loc.df$snp)]
result.df$pos<-snp.loc.df$pos[match(result.df$snps, snp.loc.df$snp)]
result.df$Symbol<-result.df$gene
nrow(result.df)
head(result.df)
max(result.df$FDR)
result.df<-result.df %>% filter(FDR<=0.05)
nrow(result.df)
output_file_name_cis
output.file<-system.file("data", package = "ProxSeek")
output.file<-paste0(output.file,"/output.csv")
output.file
if(!is.na(output.file)){
fwrite(result.df, output.file)
}else{
stop("Must have an output file!")
}
result<-ProxSeek_selectProxCis(ciseqtl.file = output.file, P_PIR_folder = P_PIR_folder, cell.type = cell.type, n.cluster = n.cluster)
#' ProxSeek_selectProxCis
#'
#' This function select cis-eQTLs that lie within promoter and promoter interaction region (P-PIR) based on tissue/cell type Hi-C data.
#' It reads cis-eQTL and cell type P-PIR files, performs filtering and annotation,
#' and returns an annotated data frame.
#'
#' @param ciseqtl.file A string specifying the path to the cis-eQTL file.
#'                     Default value is NA. It's recommended to provide the actual file path for specific analyses.
#' @param cell.type A string specifying the cell type. This is used to construct the file name of the cell type-specific P-PIR file.
#' @param n.cluster An integer specifying the number of clusters to use in parallel processing. Default is 4.
#'                  This should be set according to the hardware capabilities of the system.
#'
#' @return An annotated data frame with the results of the cis-eQTL and Hi-C interaction processing.
#'
#' @examples
#' # Example usage:
#' \dontrun{
#' demo.ciseqtl.file<-system.file("data", "ciseqtl.csv", package = "ProxSeek")
#' result <- ProxSeek_selectProxCis(ciseqtl.file = demo.ciseqtl.file, P_PIR_folder = /path/to/downloaded/PIRfiles/,
#' cell.type = "monocyte", n.cluster =4)
#' }
#'
#'
#' @import data.table
#' @import dplyr
#' @import plyr
#' @import stringr
#' @import foreach
#' @import doSNOW
#' @import parallel
#' @export
ProxSeek_selectProxCis<-function(ciseqtl.file=NA, P_PIR_folder=NA, cell.type=NA, n.cluster = 4){
###### 1. read cis-eqtl tsv file
# ciseqtl.file<-"/media/xingewang/JRlab10TB/Project_Space/Proxi_HIC_Validation/Prepare_CellTypePPIR/ciseqtl.tsv"
cis.eqtl<-fread(ciseqtl.file)
if (!file.exists(ciseqtl.file)) {
stop("Cis-eQTL file not found: ", ciseqtl.file)
}
# Attempt to read the file
tryCatch({
cis.eqtl <- fread(ciseqtl.file)
}, error = function(e) {
stop("Failed to read cis-eQTL file: ", e$message)
})
# Validate the file content
requiredColumns <- c("snps", "gene","Symbol","chr","pos")
if (!all(requiredColumns %in% colnames(cis.eqtl))) {
stop("Cis-eQTL file does not contain the required columns.")
}
###### 2. read cell.type P-PIR file
# cell.type<-"monocyte"
cell.type.file<-paste0("P_PIR_", cell.type,".csv")
cell.type.file<-paste0(P_PIR_folder,"/",cell.type.file)
# Check if the file exists
if (!file.exists(cell.type.file)) {
stop("Cell type P-PIR file not found: ", cell.type.file)
}
# Attempt to read the file
tryCatch({
cell.hic <- fread(cell.type.file)
}, error = function(e) {
stop("Failed to read cell type P-PIR file: ", e$message)
})
# Validate the file content
requiredColumns <- c("P_chr","P_start", "P_end", "P_gene","PIR_chr","PIR_start", "PIR_end", "PIR_gene")
if (!all(requiredColumns %in% colnames(cell.hic))) {
stop("Cell type P-PIR file does not contain required columns.")
}
##### 3. report genes overlapped
bait.arr<-unique(as.character(str_split(cell.hic$P_gene, ";", simplify = T)))
bait.arr<-bait.arr[bait.arr!=""]
bait.arr<-str_split(bait.arr, "[.]", simplify = T)[,1]
bait.arr<-unique(bait.arr)
hic.gene.n<-length(bait.arr)
cat(paste0("P-PIR file reported ", hic.gene.n, " genes!"))
sel.genes<-intersect(unique(cis.eqtl$Symbol), bait.arr)
share.gene.n<-length(sel.genes)
cat(paste0("Shared number of genes with cis-eQTL report ", share.gene.n))
message("")
# 4. Subset cis-eqtl with genes report
eqtl.df.sub<-cis.eqtl %>% filter(Symbol %in% sel.genes)
eqtl.df.remain<-cis.eqtl %>% filter(!Symbol %in% sel.genes)
# 5. annotate the P PIR for each SNP
##### get all the gene names into each row and split gene symbol and remove version
kk<-str_split(cell.hic$P_gene, ";", simplify = T)
kk.new<-gsub("\\..*", "", kk)
head(kk.new)
######## loop each gene and subset eQTLs
g.all<-unique(eqtl.df.sub$Symbol)
g.all<-g.all[g.all!=""]
cl <- makeCluster(n.cluster)
registerDoSNOW(cl)
iterations <- length(g.all)
pb <- txtProgressBar(max = iterations, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
tmp<-foreach(gi=g.all,.combine = 'rbind', .verbose = T, .options.snow = opts) %dopar% {
### subset eqtl for this gene
gi.eqtl<-eqtl.df.sub[eqtl.df.sub$Symbol==gi,]
### get promoter and PIR region
sel.idx<-which(rowSums(kk.new == gi)!=0)
gi.hic.df<-cell.hic[sel.idx,]
# loop each snps for this gene
in.pir<-sapply(gi.eqtl$pos, function(x){
x %in% unlist(Map(`:`, gi.hic.df$PIR_start, gi.hic.df$PIR_end))})
snp.idx.pir<-which(in.pir == 1)
in.prm <-sapply(gi.eqtl$pos, function(x){
x %in% unlist(Map(`:`, gi.hic.df$P_start, gi.hic.df$P_end))})
snp.idx.prm<-which(in.prm == 1)
#### within snp
gi.eqtl$Within<-"non_PPIR"
gi.eqtl$Within[snp.idx.pir]<-"PIR"
gi.eqtl$Within[snp.idx.prm]<-"Promoter"
return(gi.eqtl)
}
stopCluster(cl)
eqtl.df.remain$Within<-"No_HiC_Support"
cis.eqtl.annotate<-rbind(tmp, eqtl.df.remain)
return(cis.eqtl.annotate)
}
result<-ProxSeek_selectProxCis(ciseqtl.file = output.file, P_PIR_folder = P_PIR_folder, cell.type = cell.type, n.cluster = n.cluster)
demo.P.PIR<-system.file("data", package = "ProxSeek")
demo.P.PIR
P_PIR_folder<-demo.P.PIR
result<-ProxSeek_selectProxCis(ciseqtl.file = output.file, P_PIR_folder = P_PIR_folder, cell.type = cell.type, n.cluster = n.cluster)
n.cluster<-10
library(snow)
library(parallel)
result<-ProxSeek_selectProxCis(ciseqtl.file = output.file, P_PIR_folder = P_PIR_folder, cell.type = cell.type, n.cluster = n.cluster)
library(foreach)
library(doSNOW)
result<-ProxSeek_selectProxCis(ciseqtl.file = output.file, P_PIR_folder = P_PIR_folder, cell.type = cell.type, n.cluster = n.cluster)
head(result)
getwd()
devtools::document()
devtools::build()
pkgdown::build_site()
devtools::build()
devtools::document()
devtools::build_manual()
pkgdown::build_site()