MatrixQTL_implmentation_April26

matrixQTL_April26.R

require(data.table)
require(MatrixEQTL)
require(readr)
totalEM_DOS = as.matrix(totalEM_DOS[-c(1:3, 358)])
rownames(totalEM_DOS) = totalEM_DOS_SNP_POS$SNP
snps = SlicedData$new()
snps$CreateFromMatrix(totalEM_DOS)

maf.list = vector('list', length(snps))
for(sl in 1:length(snps)) {
  slice = snps[[sl]];
  maf.list[[sl]] = rowMeans(slice,na.rm=TRUE)/2;
  
  maf.list[[sl]] = pmin(maf.list[[sl]],1-maf.list[[sl]]);
}
maf = unlist(maf.list);
snps$RowReorder(maf>0.05);
show(snps)
slice = NULL
tmpf <- tempfile()
me =Matrix_eQTL_main(snps = snps, 
                     gene = EM_SERUM_PHENOS, 
                     cvrt = finalCOVS, 
                     output_file_name = tmpf,
                     pvOutputThreshold = 1e-5,
                     useModel = modelLINEAR, 
                     errorCovariance = numeric(0), 
                     verbose = T, 
                     #output_file_name.cis = output_file_name_cis,
                     #pvOutputThreshold.cis = pvOutputThreshold_cis,
                     #snpspos = chrspos, 
                     #genepos = protpos2,
                     #cisDist = cisDist,
                     pvalue.hist = "qqplot",
                     min.pv.by.genesnp = F,
                     noFDRsaveMemory = F)
unlink(tmpf)


unlink(tmpf)
total_EM_Serum_QTLS = me$all$eqtls
total_EM_Serum_QTLS = merge.data.frame(total_EM_Serum_QTLS, totalEM_DOS_SNP_POS, by.x = "snps", by.y = "SNP", all.x=T)
total_EM_Serum_QTLS= arrange(total_EM_Serum_QTLS, pvalue)
require(qqman)
manhattan(total_EM_Serum_QTLS, chr = "CHR", bp = "POS", snp = "snps", p = "pvalue")









require(data.table)
require(MatrixEQTL)
require(readr)
total_dos = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_Dosages_April24.csv', header = T, verbose = T)
write.csv(total_dos[, 1:3, with =F], file = '/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/SNP_POS_April24.csv') ### write out the snp_pos_chr so 
inversion=total_dos[CHR==8 & POS > 8135000 & POS < 11930000] ## remove chr8 pos 8.13mb to 11.93mb 
total_dos[, EM_13063 := NULL]  ### remove ataxia patient EM_13063
EM_SVDs = total_dos[!total_dos$SNP %in% inversion$SNP] %>% as.data.frame()## remove the inversion snps and calcaulate the covars
##get mafs for each snps
snps_svd = SlicedData$new()
snps_svd$CreateFromMatrix(as.matrix(EM_SVDs[-c(1:3)]))
rownames(snps_svd) = EM_SVDs$SNP
maf.list = vector('list', length(snps_svd))
for(sl in 1:length(snps_svd)) {
  slice = snps_svd[[sl]];
  maf.list[[sl]] = rowMeans(slice,na.rm=TRUE)/2;
  maf.list[[sl]] = pmin(maf.list[[sl]],1-maf.list[[sl]]);
}
maf = unlist(maf.list);
snps_svd$RowReorder(maf>0.05);
show(snps)
slice = NULL
total_dos = data.frame(row.names = rownames(snps_svd))
for (sl in 1:length(snps_svd)){
  total_dos = 
}


dupsnps=total_dos$SNP[duplicated(total_dos$SNP)]
duplicatedsnps=total_dos[total_dos$SNP %in% dupsnps]
meltdup=melt(duplicatedsnps, id.vars = c("SNP", "POS", "CHR"))
meltdup
meltdup[, list(value = mean(value)), by=list(SNP, POS, CHR, variable)]
aggmelt=meltdup[, list(value = mean(value)), by=list(SNP, POS, CHR, variable)]
aggdcast=dcast(aggmelt, formula = SNP+CHR+POS~variable, value.var = "value")
### remove duplicates and rebind the aggdcsat from above
total_dos = total_dos[!duplicated(total_dos$SNP)]




snpsrows = total_dos$SNP ## assign snps to a vector
total_dos = as.matrix(total_dos[, -c(1:3), with =F])
rownames(total_dos) = snpsrows ## assing the above vector to matrix for use as rownames

###read in the phnotype files
csf_pheno = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_CSF_Pheno_T_April24.csv', header = T)
serum_pheno = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_SERUM_Pheno_T_April24.csv', header = T)
combined_pheno = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_combined_Pheno_T_April24.csv', header = T)
ratios = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_CSF_SERUM_Ratios_T_April24.csv', header = T)
covars = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/final_COVS_April24.csv', header = T)
mafs = fread('/media/labcomp/HDD3/GWAS/final_dos/Re_analysis_April2017/EM_Maf_April24.csv', header = T)
mafs2 = mafs[maf>0.05]