### Calculate mQTLs for all probes within schizophrenia GWAS regions identified in the PGC GWAS as input for bayesian co-localistaion analysis
### Analysis is split by chromosomes

### Copyright (C) 2016  Eilis Hannon

##    This program is free software; you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##    the Free Software Foundation; either version 2 of the License, or
##    (at your option) any later version.

##    This program is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.

##   You should have received a copy of the GNU General Public License along
##    with this program; if not, write to the Free Software Foundation, Inc.,
##    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

setwd("")

library(MatrixEQTL)
base.dir = find.package("MatrixEQTL")
useModel = modelLINEAR; 

covariates_file_name = ""; 

##load covariate data
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
cvrt$fileSliceSize = 2000;      # read file in pieces of 2,000 rows
cvrt$LoadFile(covariates_file_name);



for(chr in 1:22){
	SNP_file_name = paste("", chr, ".txt", sep = "")
	snps_location_file_name = paste("", chr, ".txt", sep = "");
	expression_file_name = paste("", chr, ".txt", sep = "")

	output_file_name = paste("", chr, ".txt", sep = "")

	## threshold of results to save
	pvOutputThreshold = 1; ### need all mQTL results for each probe
	cisDist<-500000
	errorCovariance = numeric();

	 ## load SNP data
	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 pieces of 2,000 rows
	snps$LoadFile( SNP_file_name );

	## load methylation data
	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 pieces of 2,000 rows
	gene$LoadFile( expression_file_name);

	## Run the analysis
	snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
	snpspos[,1]<-rownames(snps)
	load("") ### load all illumina 450K annotation
	probeAnnot<-probeAnnot[rownames(gene), c("NAME", "CHR", "MAPINFO", "MAPINFO")]
	probeAnnot[,4]<-probeAnnot[,4]+1
	genepos<-probeAnnot

	### run eqtls
	me = Matrix_eQTL_main(
		snps = snps,
		gene = gene,
		cvrt = cvrt,
		output_file_name = output_file_name,
		pvOutputThreshold = 0,
		output_file_name.cis = output_file_name,
		pvOutputThreshold.cis = pvOutputThreshold,
		snpspos = snpspos, 
		genepos = genepos,
		cisDist = cisDist,
		useModel = useModel, 
		errorCovariance = errorCovariance, 
		verbose = TRUE,
		pvalue.hist = "qqplot",
		min.pv.by.genesnp = FALSE,
		noFDRsaveMemory = FALSE)
		
}