rm(list = ls())
gc()
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 372365 19.9 750400 40.1 592000 31.7
## Vcells 629976 4.9 1308461 10.0 937336 7.2
# set the working directory
setwd("/Volumes/Transcend/Thesis_project/subsetted_liver")
# subset dataset
sebsetn <- 30
mouse.liver.expression.eqtl <-read.table(file="2016-05-16 mouse.liver.expression.eqtl.txt", header=T)
head(mouse.liver.expression.eqtl)
## ProbeSet BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18
## 1 1415670_at 10.090 10.200 10.300 10.208 9.830 10.238 9.914 10.348
## 2 1415671_at 10.932 11.088 11.007 11.020 10.955 11.120 11.012 11.123
## 3 1415672_at 11.432 11.417 11.442 11.555 11.561 11.318 11.461 11.561
## 4 1415673_at 7.535 7.382 7.566 7.162 7.403 7.342 7.213 7.581
## 5 1415674_a_at 9.757 9.972 9.269 9.873 9.354 9.918 9.459 9.655
## 6 1415675_at 9.029 9.009 9.245 9.282 9.415 9.098 9.060 8.937
## BXD19 BXD2 BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31
## 1 9.939 9.871 10.077 10.159 9.746 9.890 10.286 10.177 9.959 9.882
## 2 10.922 10.802 10.988 10.969 11.104 10.979 10.905 11.013 11.071 10.996
## 3 11.575 11.426 11.367 11.328 11.499 11.531 11.666 11.500 11.502 11.446
## 4 7.551 7.368 7.251 7.373 7.408 7.473 7.238 7.424 7.475 7.491
## 5 9.544 9.557 9.460 9.322 9.795 9.758 9.628 9.263 9.671 9.645
## 6 8.995 9.104 9.123 9.086 9.087 9.040 8.844 9.173 8.951 9.131
## BXD32 BXD33 BXD34 BXD36 BXD38 BXD39 BXD40 BXD42 BXD5 BXD6
## 1 10.059 10.102 10.174 10.022 10.364 9.745 10.074 9.961 10.160 10.069
## 2 10.854 11.084 11.059 10.923 11.053 11.030 11.067 10.891 10.878 11.043
## 3 11.548 11.511 11.490 11.545 11.496 11.516 11.543 11.457 11.374 11.504
## 4 7.419 7.386 7.496 7.369 7.265 7.130 7.149 7.383 7.439 7.233
## 5 9.368 9.741 9.721 9.684 9.387 9.649 9.693 9.793 9.840 9.519
## 6 9.108 9.122 9.153 9.008 9.138 9.126 9.233 9.060 9.134 9.081
## BXD8 BXD9
## 1 9.956 10.142
## 2 11.206 10.982
## 3 11.531 11.700
## 4 7.397 7.390
## 5 9.340 9.160
## 6 8.926 9.274
dim(mouse.liver.expression.eqtl)
## [1] 20634 31
set.seed(50)
sub.mouse.liver.expression.eqtl <- mouse.liver.expression.eqtl[, c(1, sample(2:dim(mouse.liver.expression.eqtl)[2],sebsetn, replace=FALSE))]
head(sub.mouse.liver.expression.eqtl)
## ProbeSet BXD36 BXD24 BXD15 BXD34 BXD24a BXD11 BXD29 BXD27
## 1 1415670_at 10.022 9.746 10.238 10.174 9.890 10.200 9.959 10.286
## 2 1415671_at 10.923 11.104 11.120 11.059 10.979 11.088 11.071 10.905
## 3 1415672_at 11.545 11.499 11.318 11.490 11.531 11.417 11.502 11.666
## 4 1415673_at 7.369 7.408 7.342 7.496 7.473 7.382 7.475 7.238
## 5 1415674_a_at 9.684 9.795 9.918 9.721 9.758 9.972 9.671 9.628
## 6 1415675_at 9.008 9.087 9.098 9.153 9.040 9.009 8.951 8.844
## BXD1 BXD12 BXD18 BXD6 BXD21 BXD40 BXD14 BXD20 BXD31 BXD28
## 1 10.090 10.300 10.348 10.069 10.159 10.074 9.830 10.077 9.882 10.177
## 2 10.932 11.007 11.123 11.043 10.969 11.067 10.955 10.988 10.996 11.013
## 3 11.432 11.442 11.561 11.504 11.328 11.543 11.561 11.367 11.446 11.500
## 4 7.535 7.566 7.581 7.233 7.373 7.149 7.403 7.251 7.491 7.424
## 5 9.757 9.269 9.655 9.519 9.322 9.693 9.354 9.460 9.645 9.263
## 6 9.029 9.245 8.937 9.081 9.086 9.233 9.415 9.123 9.131 9.173
## BXD9 BXD39 BXD16 BXD38 BXD32 BXD33 BXD19 BXD13 BXD5 BXD42
## 1 10.142 9.745 9.914 10.364 10.059 10.102 9.939 10.208 10.160 9.961
## 2 10.982 11.030 11.012 11.053 10.854 11.084 10.922 11.020 10.878 10.891
## 3 11.700 11.516 11.461 11.496 11.548 11.511 11.575 11.555 11.374 11.457
## 4 7.390 7.130 7.213 7.265 7.419 7.386 7.551 7.162 7.439 7.383
## 5 9.160 9.649 9.459 9.387 9.368 9.741 9.544 9.873 9.840 9.793
## 6 9.274 9.126 9.060 9.138 9.108 9.122 8.995 9.282 9.134 9.060
## BXD8 BXD2
## 1 9.956 9.871
## 2 11.206 10.802
## 3 11.531 11.426
## 4 7.397 7.368
## 5 9.340 9.557
## 6 8.926 9.104
dim(sub.mouse.liver.expression.eqtl)
## [1] 20634 31
write.table(sub.mouse.liver.expression.eqtl,file="2016-05-09 sub.mouse.liver.expression.eqtl.txt", sep="\t", row.names=FALSE, quote=FALSE)
#subset liver snp expression data
BXD.geno.SNP.eqtl.for.liver <-read.table(file="2016-05-16 BXD.geno.SNP.eqtl.for.liver.txt", header=T)
head(BXD.geno.SNP.eqtl.for.liver)
## Locus BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD2
## 1 rs13459051 0 1 0 1 0 1 0 0 0 0
## 2 rs13459060 0 0 0 1 1 0 1 0 1 0
## 3 rs13459062 1 0 0 1 1 0 1 0 1 0
## 4 rs13459063 0 0 0 1 1 0 1 0 1 0
## 5 rs13459064 1 1 0 1 1 0 1 1 1 1
## 6 rs13459069 0 1 1 0 0 1 0 0 0 0
## BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31 BXD32 BXD33 BXD34 BXD36
## 1 1 0 1 1 1 1 1 0 1 1 1 0
## 2 0 1 1 1 0 0 0 1 1 0 0 1
## 3 0 1 1 1 0 0 0 1 1 0 0 1
## 4 0 1 1 1 0 0 0 1 1 0 0 1
## 5 1 1 1 1 0 0 1 1 1 0 0 1
## 6 0 1 0 0 1 1 0 0 0 0 1 1
## BXD38 BXD39 BXD40 BXD42 BXD5 BXD6 BXD8 BXD9
## 1 0 1 0 0 1 1 0 0
## 2 0 0 1 1 1 0 0 1
## 3 0 0 1 1 1 0 1 1
## 4 0 0 1 1 1 0 1 1
## 5 1 1 1 1 1 0 0 0
## 6 1 1 0 1 1 0 0 1
dim(BXD.geno.SNP.eqtl.for.liver)
## [1] 3024 31
set.seed(50)
sub.BXD.geno.SNP.eqtl.for.liver <- BXD.geno.SNP.eqtl.for.liver[, c(1, sample(2:dim(BXD.geno.SNP.eqtl.for.liver)[2],sebsetn, replace=FALSE))]
head(sub.BXD.geno.SNP.eqtl.for.liver)
## Locus BXD36 BXD24 BXD15 BXD34 BXD24a BXD11 BXD29 BXD27 BXD1 BXD12
## 1 rs13459051 0 1 1 1 1 1 1 1 0 0
## 2 rs13459060 1 1 0 0 1 0 0 0 0 0
## 3 rs13459062 1 1 0 0 1 0 0 0 1 0
## 4 rs13459063 1 1 0 0 1 0 0 0 0 0
## 5 rs13459064 1 1 0 0 1 1 1 0 1 0
## 6 rs13459069 1 0 1 1 0 1 0 1 0 1
## BXD18 BXD6 BXD21 BXD40 BXD14 BXD20 BXD31 BXD28 BXD9 BXD39 BXD16 BXD38
## 1 0 1 0 0 0 1 0 1 0 1 0 0
## 2 0 0 1 1 1 0 1 0 1 0 1 0
## 3 0 0 1 1 1 0 1 0 1 0 1 0
## 4 0 0 1 1 1 0 1 0 1 0 1 0
## 5 1 0 1 1 1 1 1 0 0 1 1 1
## 6 0 0 1 0 0 0 0 1 1 1 0 1
## BXD32 BXD33 BXD19 BXD13 BXD5 BXD42 BXD8 BXD2
## 1 1 1 0 1 1 0 0 0
## 2 1 0 1 1 1 1 0 0
## 3 1 0 1 1 1 1 1 0
## 4 1 0 1 1 1 1 1 0
## 5 1 0 1 1 1 1 0 1
## 6 0 0 0 0 1 1 0 0
dim(sub.BXD.geno.SNP.eqtl.for.liver)
## [1] 3024 31
write.table(sub.BXD.geno.SNP.eqtl.for.liver,file="2016-05-09 sub.BXD.geno.SNP.eqtl.for.liver.txt", sep="\t", row.names=FALSE, quote=FALSE)
library(MatrixEQTL)
## Location of the package with the data files.
base.dir = "/Volumes/Transcend/Thesis_project/subsetted_liver";
## Settings
# Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
# Genotype file name
SNP_file_name = paste(base.dir, "/2016-05-09 sub.BXD.geno.SNP.eqtl.for.liver.txt", sep="");
snps_location_file_name = paste(base.dir, "/2016-05-16 BXD.geno.loc.eqtl.for.liver.txt", sep="");
# Gene expression file name
expression_file_name = paste(base.dir, "/2016-05-09 sub.mouse.liver.expression.eqtl.txt", sep="");
gene_location_file_name = paste(base.dir, "/2016-05-16 liver.gene.loc.txt", sep="");
# Covariates file name
# Set to character() for no covariates
covariates_file_name = character() ;
# Output file name
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();
# Only associations significant at this level will be saved
pvOutputThreshold_cis = 1;
pvOutputThreshold_tra = 0.000000000000005;
# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric();
# errorCovariance = read.table("Sample_Data/errorCovariance.txt");
# Distance for local gene-SNP pairs
cisDist = 1e6;
## Load genotype data
snps = SlicedData$new();
snps$fileDelimiter = "\t"; # the TAB character
snps$fileOmitCharacters = "NA"; # denote missing values;
snps$fileSkipRows = 1;
snps$fileSkipColumns = 1;
snps$fileSliceSize = 2000;
snps$LoadFile(SNP_file_name);
## Rows read: 2,000
## Rows read: 3024 done.
## Load gene expression data
gene = SlicedData$new();
gene$fileDelimiter = "\t";
gene$fileOmitCharacters = "NA"; # denote missing values;
gene$fileSkipRows = 1;
gene$fileSkipColumns = 1;
gene$fileSliceSize = 2000;
gene$LoadFile(expression_file_name);
## Rows read: 2,000
## Rows read: 4,000
## Rows read: 6,000
## Rows read: 8,000
## Rows read: 10,000
## Rows read: 12,000
## Rows read: 14,000
## Rows read: 16,000
## Rows read: 18,000
## Rows read: 20,000
## Rows read: 20634 done.
## Load covariates
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);
head(genepos)
## probeset Chromosome Start End
## 1 1415670_at 6 87887814 87913595
## 2 1415671_at 8 105524465 105566047
## 3 1415672_at 8 23241353 23257074
## 4 1415673_at 5 129765558 129787449
## 5 1415674_a_at 9 44403761 44407548
## 6 1415675_at 2 32570858 32573579
me = Matrix_eQTL_main(
snps = snps,
gene = gene,
output_file_name = output_file_name_tra,
pvOutputThreshold = pvOutputThreshold_tra,
useModel = useModel,
errorCovariance = numeric(),
verbose = TRUE,
output_file_name.cis = output_file_name_cis,
pvOutputThreshold.cis = pvOutputThreshold_cis,
snpspos = snpspos,
genepos = genepos,
cisDist = cisDist,
pvalue.hist = TRUE,
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
## Matching data files and location files
## 20634 of 20634 genes matched
## 3024 of 3024 SNPs matched
## Task finished in 0.017 seconds
## Reordering SNPs
##
## Task finished in 0.088 seconds
## Reordering genes
##
## Task finished in 0.142 seconds
## Processing covariates
## Task finished in 0.001 seconds
## Processing gene expression data (imputation, residualization, etc.)
## Task finished in 0.028 seconds
## Creating output file(s)
## Task finished in 0.011 seconds
## Performing eQTL analysis
## 4.54% done, 5,899 cis-eQTLs, 30 trans-eQTLs
## 9.09% done, 11,411 cis-eQTLs, 37 trans-eQTLs
## 13.63% done, 17,998 cis-eQTLs, 83 trans-eQTLs
## 18.18% done, 22,514 cis-eQTLs, 90 trans-eQTLs
## 22.72% done, 27,349 cis-eQTLs, 177 trans-eQTLs
## 27.27% done, 31,978 cis-eQTLs, 178 trans-eQTLs
## 31.81% done, 38,098 cis-eQTLs, 196 trans-eQTLs
## 36.36% done, 38,605 cis-eQTLs, 196 trans-eQTLs
## 40.90% done, 196 trans-eQTLs
## 45.45% done, 196 trans-eQTLs
## 50.00% done, 198 trans-eQTLs
## 54.54% done, 198 trans-eQTLs
## 59.09% done, 198 trans-eQTLs
## 63.63% done, 198 trans-eQTLs
## 68.18% done, 198 trans-eQTLs
## 72.72% done, 206 trans-eQTLs
## 77.27% done, 206 trans-eQTLs
## 81.81% done, 206 trans-eQTLs
## 86.36% done, 42,998 cis-eQTLs, 222 trans-eQTLs
## 90.90% done, 48,565 cis-eQTLs, 232 trans-eQTLs
## 95.45% done, 54,309 cis-eQTLs, 258 trans-eQTLs
## 100.00% done, 54,955 cis-eQTLs, 263 trans-eQTLs
## Task finished in 9.549 seconds
##
unlink(output_file_name_cis);
## Results:
cat('Analysis done in:', me$time.in.sec, ' seconds', '\n')
## Analysis done in: 9.109 seconds
cat('Detected local eQTLs:','\n')
## Detected local eQTLs:
cis.eqtls<-me$cis$eqtls
head(cis.eqtls)
## snps gene statistic pvalue FDR beta
## 1 rs4163042 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
## 2 rs4163058 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
## 3 rs4163391 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
## 4 rs4151923 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
## 5 rs3090019 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
## 6 rs4164131 1452705_at -76.34025 4.891311e-34 3.840028e-30 -4.190552
dim(cis.eqtls)
## [1] 54955 6
cis.eqtls$beta_se <-cis.eqtls$beta/cis.eqtls$statistic
write.table(cis.eqtls,file="2016-05-09 sub.mouseliver.cis.1M.eqtls.txt", sep="\t", row.names=FALSE, quote=FALSE)
# eqtl analysis for lung
## Location of the package with the data files.
base.dir = "/Volumes/Transcend/Thesis_project/subsetted_liver";
## Settings
# Linear model to use, modelANOVA, modelLINEAR, or modelLINEAR_CROSS
useModel = modelLINEAR; # modelANOVA, modelLINEAR, or modelLINEAR_CROSS
# Genotype file name
SNP_file_name = paste(base.dir, "/2016-05-16 BXD.geno.SNP.eqtl.for.lung.txt", sep="");
snps_location_file_name = paste(base.dir, "/2016-05-16 BXD.geno.loc.eqtl.for.lung.txt", sep="");
# Gene expression file name
expression_file_name = paste(base.dir, "/2016-05-16 mouse.lung.expression.eqtl.txt", sep="");
gene_location_file_name = paste(base.dir, "/2016-05-16 lung.gene.loc.txt", sep="");
# Covariates file name
# Set to character() for no covariates
covariates_file_name = character() ;
# Output file name
output_file_name_cis = tempfile();
output_file_name_tra = tempfile();
# Only associations significant at this level will be saved
pvOutputThreshold_cis = 1;
pvOutputThreshold_tra = 0.000000000000005;
# Error covariance matrix
# Set to numeric() for identity.
errorCovariance = numeric();
# errorCovariance = read.table("Sample_Data/errorCovariance.txt");
# Distance for local gene-SNP pairs
cisDist = 1e6;
## Load genotype data
snps = SlicedData$new();
snps$fileDelimiter = "\t"; # the TAB character
snps$fileOmitCharacters = "NA"; # denote missing values;
snps$fileSkipRows = 1;
snps$fileSkipColumns = 1;
snps$fileSliceSize = 2000;
snps$LoadFile(SNP_file_name);
## Rows read: 2,000
## Rows read: 3024 done.
## Load gene expression data
gene = SlicedData$new();
gene$fileDelimiter = "\t";
gene$fileOmitCharacters = "NA"; # denote missing values;
gene$fileSkipRows = 1;
gene$fileSkipColumns = 1;
gene$fileSliceSize = 2000;
gene$LoadFile(expression_file_name);
## Rows read: 2,000
## Rows read: 4,000
## Rows read: 6,000
## Rows read: 8,000
## Rows read: 10,000
## Rows read: 12,000
## Rows read: 14,000
## Rows read: 16,000
## Rows read: 18,000
## Rows read: 20,000
## Rows read: 20634 done.
## Load covariates
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);
head(genepos)
## probeset Chromosome Start End
## 1 1415670_at 6 87887814 87913595
## 2 1415671_at 8 105524465 105566047
## 3 1415672_at 8 23241353 23257074
## 4 1415673_at 5 129765558 129787449
## 5 1415674_a_at 9 44403761 44407548
## 6 1415675_at 2 32570858 32573579
me = Matrix_eQTL_main(
snps = snps,
gene = gene,
output_file_name = output_file_name_tra,
pvOutputThreshold = pvOutputThreshold_tra,
useModel = useModel,
errorCovariance = numeric(),
verbose = TRUE,
output_file_name.cis = output_file_name_cis,
pvOutputThreshold.cis = pvOutputThreshold_cis,
snpspos = snpspos,
genepos = genepos,
cisDist = cisDist,
pvalue.hist = TRUE,
min.pv.by.genesnp = FALSE,
noFDRsaveMemory = FALSE);
## Matching data files and location files
## 20634 of 20634 genes matched
## 3024 of 3024 SNPs matched
## Task finished in 0.026 seconds
## Reordering SNPs
##
## Task finished in 0.352 seconds
## Reordering genes
##
## Task finished in 0.336 seconds
## Processing covariates
## Task finished in 0.002 seconds
## Processing gene expression data (imputation, residualization, etc.)
## Task finished in 0.06 seconds
## Creating output file(s)
## Task finished in 0.012 seconds
## Performing eQTL analysis
## 4.54% done, 5,899 cis-eQTLs, 87 trans-eQTLs
## 9.09% done, 11,411 cis-eQTLs, 231 trans-eQTLs
## 13.63% done, 17,998 cis-eQTLs, 287 trans-eQTLs
## 18.18% done, 22,514 cis-eQTLs, 326 trans-eQTLs
## 22.72% done, 27,349 cis-eQTLs, 422 trans-eQTLs
## 27.27% done, 31,978 cis-eQTLs, 452 trans-eQTLs
## 31.81% done, 38,098 cis-eQTLs, 511 trans-eQTLs
## 36.36% done, 38,605 cis-eQTLs, 514 trans-eQTLs
## 40.90% done, 514 trans-eQTLs
## 45.45% done, 514 trans-eQTLs
## 50.00% done, 516 trans-eQTLs
## 54.54% done, 518 trans-eQTLs
## 59.09% done, 518 trans-eQTLs
## 63.63% done, 518 trans-eQTLs
## 68.18% done, 518 trans-eQTLs
## 72.72% done, 518 trans-eQTLs
## 77.27% done, 518 trans-eQTLs
## 81.81% done, 518 trans-eQTLs
## 86.36% done, 42,998 cis-eQTLs, 577 trans-eQTLs
## 90.90% done, 48,565 cis-eQTLs, 615 trans-eQTLs
## 95.45% done, 54,309 cis-eQTLs, 663 trans-eQTLs
## 100.00% done, 54,955 cis-eQTLs, 671 trans-eQTLs
## Task finished in 9.303 seconds
##
unlink(output_file_name_cis);
## Results:
cat('Analysis done in:', me$time.in.sec, ' seconds', '\n')
## Analysis done in: 9.473 seconds
cat('Detected local eQTLs:','\n')
## Detected local eQTLs:
cis.eqtls<-me$cis$eqtls
head(cis.eqtls)
## snps gene statistic pvalue FDR beta
## 1 rs3681670 1421144_at -75.53535 2.361797e-47 6.489628e-43 -3.093200
## 2 rs3666583 1421144_at -75.53535 2.361797e-47 6.489628e-43 -3.093200
## 3 rs13481364 1438758_at -50.66421 5.603003e-40 4.398757e-36 -2.167143
## 4 rs3708704 1438758_at -50.66421 5.603003e-40 4.398757e-36 -2.167143
## 5 rs13481365 1438758_at -50.66421 5.603003e-40 4.398757e-36 -2.167143
## 6 rs3685299 1438758_at -50.66421 5.603003e-40 4.398757e-36 -2.167143
dim(cis.eqtls)
## [1] 54955 6
cis.eqtls$beta_se <-cis.eqtls$beta/cis.eqtls$statistic
write.table(cis.eqtls,file="2016-05-16 mouselung.cis.1M.eqtls.txt", sep="\t", row.names=FALSE, quote=FALSE)
# load mouse lung cis eqtl result
lung.mouse.eQTL<-read.table(file="2016-05-16 mouselung.cis.1M.eqtls.txt", header=T)
# load mouse liver cis eqtl result
liver.mouse.eQTL<-read.table(file="2016-05-09 sub.mouseliver.cis.1M.eqtls.txt", header=T)
mouse4302ensembl_id<-read.table(file="2015-12-04 mouse4302ensembl_id.txt", header=T)
mouse430aensembl_id<-read.table(file="2015-12-07 mouse430aensembl_id.txt", header=T)
# Add ensemble id annoatation to the data
lung.mouse.eQTL<-merge(lung.mouse.eQTL, mouse4302ensembl_id, by.x = "gene", by.y="probe_id")
liver.mouse.eQTL<-merge(liver.mouse.eQTL, mouse430aensembl_id, by.x = "gene", by.y="probe_id")
head(lung.mouse.eQTL)
## gene snps statistic pvalue FDR beta
## 1 1415670_at rs13478880 1.1363521 0.2621028 0.5270176 0.06786667
## 2 1415670_at rs13478876 0.9133601 0.3661462 0.6302137 0.05583410
## 3 1415670_at rs3713705 1.1363521 0.2621028 0.5270176 0.06786667
## 4 1415670_at rs13475374 1.0503967 0.2994035 0.5651655 0.06286667
## 5 1415672_at rs13479651 -1.3052534 0.1987480 0.4483475 -0.03993582
## 6 1415672_at rs3661760 -1.3052534 0.1987480 0.4483475 -0.03993582
## beta_se ensembl_id
## 1 0.05972327 ENSMUSG00000030058
## 2 0.06113044 ENSMUSG00000030058
## 3 0.05972327 ENSMUSG00000030058
## 4 0.05985040 ENSMUSG00000030058
## 5 0.03059622 ENSMUSG00000015341
## 6 0.03059622 ENSMUSG00000015341
head(liver.mouse.eQTL)
## gene snps statistic pvalue FDR beta
## 1 1415670_at rs13475374 0.1359233 0.8928544 0.9656043 0.008533333
## 2 1415670_at rs13478880 0.3577236 0.7232327 0.9020022 0.022464286
## 3 1415670_at rs13478876 0.3577236 0.7232327 0.9020022 0.022464286
## 4 1415670_at rs3713705 0.3577236 0.7232327 0.9020022 0.022464286
## 5 1415672_at rs13479651 -0.1440997 0.8864539 0.9638256 -0.004607143
## 6 1415672_at rs3661760 -0.1440997 0.8864539 0.9638256 -0.004607143
## beta_se ensembl_id
## 1 0.06278048 ENSMUSG00000030058
## 2 0.06279788 ENSMUSG00000030058
## 3 0.06279788 ENSMUSG00000030058
## 4 0.06279788 ENSMUSG00000030058
## 5 0.03197192 ENSMUSG00000015341
## 6 0.03197192 ENSMUSG00000015341
library(data.table)
# Select Gene-SNP pair with minimum P value
lung.mouse.eQTL.min <- data.table(lung.mouse.eQTL, key=c('ensembl_id', "pvalue"))
lung.mouse.eQTL.min<-lung.mouse.eQTL.min[J(unique(ensembl_id)),mult="first"]
lung.mouse.eQTL.min<-as.data.frame(lung.mouse.eQTL.min)
liver.mouse.eQTL.min <- data.table(liver.mouse.eQTL, key=c('ensembl_id', "pvalue"))
liver.mouse.eQTL.min<-liver.mouse.eQTL.min[J(unique(ensembl_id)),mult="first"]
liver.mouse.eQTL.min<-as.data.frame(liver.mouse.eQTL.min)
library(plyr)
lung.mouse.eQTL.min<-rename(lung.mouse.eQTL.min, c("pvalue"="lung_pvalue", "beta"="lung.beta", "beta_se"="lung.beta_se"))
liver.mouse.eQTL.min<-rename(liver.mouse.eQTL.min, c("pvalue"="liver_pvalue", "beta"="liver.beta", "beta_se"="liver.beta_se"))
head(lung.mouse.eQTL.min)
## gene snps statistic lung_pvalue FDR lung.beta
## 1 1428645_at rs6259798 -1.4324741 0.159238340 0.39180403 -0.05980000
## 2 1416677_at rs3707743 -0.1308388 0.896513064 0.95927760 -0.02366071
## 3 1451678_at rs13481264 -0.3904599 0.698124238 0.86317736 -0.01356126
## 4 1425955_at rs13478643 -1.0185410 0.314117648 0.58046119 -0.07813675
## 5 1426241_a_at rs4224744 0.9192838 0.363075471 0.62709198 0.05284706
## 6 1439267_x_at rs4136521 3.4048233 0.001444641 0.01001896 0.18997619
## lung.beta_se ensembl_id
## 1 0.04174596 ENSMUSG00000000001
## 2 0.18083869 ENSMUSG00000000049
## 3 0.03473152 ENSMUSG00000000056
## 4 0.07671439 ENSMUSG00000000058
## 5 0.05748721 ENSMUSG00000000085
## 6 0.05579620 ENSMUSG00000000088
head(liver.mouse.eQTL.min)
## gene snps statistic liver_pvalue FDR liver.beta
## 1 1428645_at rs6259798 -0.6238563 5.377718e-01 0.8121255718 -0.03197222
## 2 1416677_at rs3704454 5.4086629 9.084753e-06 0.0001768518 0.08446667
## 3 1451677_at rs3726373 1.0590703 2.986190e-01 0.6348151014 0.06783732
## 4 1425955_at rs3023064 -1.9727688 5.846922e-02 0.2482175362 -0.05122222
## 5 1426241_a_at rs3675629 -2.1399223 4.121599e-02 0.1944895183 -0.06860287
## 6 1448153_at rs13480217 5.1217929 1.987922e-05 0.0003521801 0.19240670
## liver.beta_se ensembl_id
## 1 0.05124934 ENSMUSG00000000001
## 2 0.01561692 ENSMUSG00000000049
## 3 0.06405365 ENSMUSG00000000056
## 4 0.02596463 ENSMUSG00000000058
## 5 0.03205858 ENSMUSG00000000085
## 6 0.03756628 ENSMUSG00000000088
tail(liver.mouse.eQTL.min)
## gene snps statistic liver_pvalue FDR
## 10596 1434694_at rs13459062 -2.3077250 0.02861902 0.1506185
## 10597 1437645_at rs13482752 2.0047987 0.05474021 0.2367767
## 10598 1449939_s_at rs13481643 0.6744392 0.50556296 0.7933074
## 10599 1422547_at rs4165065 -1.6676237 0.10653819 0.3616534
## 10600 1451476_at rs4165065 0.3972613 0.69418700 0.8906254
## 10601 1448115_at rs4165069 1.5491597 0.13257346 0.4121266
## liver.beta liver.beta_se ensembl_id
## 10596 -0.12408929 0.05377126 ENSMUSG00000099041
## 10597 0.06022624 0.03004104 ENSMUSG00000099083
## 10598 0.01886111 0.02796562 ENSMUSG00000099116
## 10599 -0.07874163 0.04721786 ENSMUSG00000099164
## 10600 0.01569378 0.03950493 ENSMUSG00000099262
## 10601 0.08752778 0.05650017 ENSMUSG00000099305
dim(lung.mouse.eQTL.min)
## [1] 10601 8
dim(liver.mouse.eQTL.min)
## [1] 10601 8
# lung, liver eqtl with ensemble_id
merged.mouse.eQTL.min<-merge(lung.mouse.eQTL.min, liver.mouse.eQTL.min, by.x = "ensembl_id", by.y="ensembl_id")
head(merged.mouse.eQTL.min)
## ensembl_id gene.x snps.x statistic.x lung_pvalue
## 1 ENSMUSG00000000001 1428645_at rs6259798 -1.4324741 0.159238340
## 2 ENSMUSG00000000049 1416677_at rs3707743 -0.1308388 0.896513064
## 3 ENSMUSG00000000056 1451678_at rs13481264 -0.3904599 0.698124238
## 4 ENSMUSG00000000058 1425955_at rs13478643 -1.0185410 0.314117648
## 5 ENSMUSG00000000085 1426241_a_at rs4224744 0.9192838 0.363075471
## 6 ENSMUSG00000000088 1439267_x_at rs4136521 3.4048233 0.001444641
## FDR.x lung.beta lung.beta_se gene.y snps.y statistic.y
## 1 0.39180403 -0.05980000 0.04174596 1428645_at rs6259798 -0.6238563
## 2 0.95927760 -0.02366071 0.18083869 1416677_at rs3704454 5.4086629
## 3 0.86317736 -0.01356126 0.03473152 1451677_at rs3726373 1.0590703
## 4 0.58046119 -0.07813675 0.07671439 1425955_at rs3023064 -1.9727688
## 5 0.62709198 0.05284706 0.05748721 1426241_a_at rs3675629 -2.1399223
## 6 0.01001896 0.18997619 0.05579620 1448153_at rs13480217 5.1217929
## liver_pvalue FDR.y liver.beta liver.beta_se
## 1 5.377718e-01 0.8121255718 -0.03197222 0.05124934
## 2 9.084753e-06 0.0001768518 0.08446667 0.01561692
## 3 2.986190e-01 0.6348151014 0.06783732 0.06405365
## 4 5.846922e-02 0.2482175362 -0.05122222 0.02596463
## 5 4.121599e-02 0.1944895183 -0.06860287 0.03205858
## 6 1.987922e-05 0.0003521801 0.19240670 0.03756628
dim(merged.mouse.eQTL.min)
## [1] 10601 15
merged.mouse.eQTL.min<-data.frame(merged.mouse.eQTL.min)
merged.mouse.eQTL.min<-merged.mouse.eQTL.min[, c(1, 5, 7, 8, 12, 14, 15 )]
head(merged.mouse.eQTL.min)
## ensembl_id lung_pvalue lung.beta lung.beta_se liver_pvalue
## 1 ENSMUSG00000000001 0.159238340 -0.05980000 0.04174596 5.377718e-01
## 2 ENSMUSG00000000049 0.896513064 -0.02366071 0.18083869 9.084753e-06
## 3 ENSMUSG00000000056 0.698124238 -0.01356126 0.03473152 2.986190e-01
## 4 ENSMUSG00000000058 0.314117648 -0.07813675 0.07671439 5.846922e-02
## 5 ENSMUSG00000000085 0.363075471 0.05284706 0.05748721 4.121599e-02
## 6 ENSMUSG00000000088 0.001444641 0.18997619 0.05579620 1.987922e-05
## liver.beta liver.beta_se
## 1 -0.03197222 0.05124934
## 2 0.08446667 0.01561692
## 3 0.06783732 0.06405365
## 4 -0.05122222 0.02596463
## 5 -0.06860287 0.03205858
## 6 0.19240670 0.03756628
write.table(merged.mouse.eQTL.min,file="2016-05-09 mouse.liver.expression.min.txt", sep="\t", row.names=FALSE, quote=FALSE)
merged.mouse.eQTL.min.variance2<-read.table(file="2016-05-09 mouse.liver.expression.min.txt", header=T)
head(merged.mouse.eQTL.min.variance2)
## ensembl_id lung_pvalue lung.beta lung.beta_se liver_pvalue
## 1 ENSMUSG00000000001 0.159238340 -0.05980000 0.04174596 5.377718e-01
## 2 ENSMUSG00000000049 0.896513064 -0.02366071 0.18083869 9.084753e-06
## 3 ENSMUSG00000000056 0.698124238 -0.01356126 0.03473152 2.986190e-01
## 4 ENSMUSG00000000058 0.314117648 -0.07813675 0.07671439 5.846922e-02
## 5 ENSMUSG00000000085 0.363075471 0.05284706 0.05748721 4.121599e-02
## 6 ENSMUSG00000000088 0.001444641 0.18997619 0.05579620 1.987922e-05
## liver.beta liver.beta_se
## 1 -0.03197222 0.05124934
## 2 0.08446667 0.01561692
## 3 0.06783732 0.06405365
## 4 -0.05122222 0.02596463
## 5 -0.06860287 0.03205858
## 6 0.19240670 0.03756628
# caculate the absolute value of live/lung beta
merged.mouse.eQTL.min.variance2$abs_liver.beta<-abs(merged.mouse.eQTL.min.variance2$liver.beta)
merged.mouse.eQTL.min.variance2$abs_lung.beta<-abs(merged.mouse.eQTL.min.variance2$lung.beta)
# caculate negative log lung p value
merged.mouse.eQTL.min.variance2$neg_log_lung_pvalue<--log10(merged.mouse.eQTL.min.variance2$lung_pvalue)
# Simple linear regression between abs_liver.beta and abs_lung.beta
# fit1<-summary(lm(abs_liver.beta ~ abs_lung.beta, data=merged.mouse.eQTL.min.variance2))
# fit1
# tau<-fit1$sigma**2
# check association between abs_liver.beta and abs.lung.beta
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
ggplot(merged.mouse.eQTL.min.variance2, aes(x=abs_lung.beta, y=abs_liver.beta)) +geom_point()+geom_smooth(method=lm)
cor(merged.mouse.eQTL.min.variance2$abs_lung.beta, merged.mouse.eQTL.min.variance2$abs_liver.beta)
## [1] 0.3644936
merged.mouse.eQTL<-merged.mouse.eQTL.min.variance2
# retrieve ensembl_id
markers<-merged.mouse.eQTL[, 1]
# Yg=Ag + Bg*Xsnp+V
# retrieve betas.hat (liver.beta)
betas.hat<-merged.mouse.eQTL$abs_liver.beta
# retrieve liver.beta_se
se<-merged.mouse.eQTL$liver.beta_se
# creat Z matrix with 2 columns: 1 for intercept,abs_lung.beta (merged.mouse.eQTL[,10])
Z<-as.matrix(merged.mouse.eQTL$abs_lung.beta)
Z<-replace(Z,is.na(Z),0)
Z<-data.frame(1,Z)
Z<-as.matrix(Z)
rowLength<-length(markers)
# liver.betas=Z*gama+T^2
# Regression: abs_liver.beta = intercept + beta*abs_lung.beta + error
lmsummary<-summary(lm(abs_liver.beta~-1+Z, data=merged.mouse.eQTL))
lmsummary
##
## Call:
## lm(formula = abs_liver.beta ~ -1 + Z, data = merged.mouse.eQTL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1632 -0.0660 -0.0403 0.0102 4.4094
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## ZX1 0.062709 0.002242 27.98 <2e-16 ***
## ZZ 0.346325 0.008594 40.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1977 on 10599 degrees of freedom
## Multiple R-squared: 0.3145, Adjusted R-squared: 0.3144
## F-statistic: 2431 on 2 and 10599 DF, p-value: < 2.2e-16
# error ~ N(0, Tau)
tau<-lmsummary$sigma**2
tau
## [1] 0.03909984
# output coeffieients (gamma matrix)
# gamma matrix
gamma<-as.matrix(lmsummary$coefficients[,1])
# trasnpose Z matrix
Z_transpose<-t(Z)
# create identity matrix
identity<-diag(nrow=rowLength)
# original betas.hat
betas.hat<-as.matrix(betas.hat)
#creat V matrix for liver_residual_variance
V <- matrix(0, rowLength, rowLength)
# V, liver residual variance
diag(V) <- merged.mouse.eQTL$liver.beta_se^2
# Creat Tau matrix
Tau<- diag(tau, rowLength, rowLength)
# follow Chen's paper and cacualte s
s <-V + Tau
# create inverse function for inversing diagnoal matrix
diag.inverse <- function(x){diag(1/diag(x), nrow(x), ncol(x))}
# create multiplication function for multiplicating two diagnoal matrix
diag.multi <- function(x,y){diag(diag(x)*diag(y), nrow(x), ncol(x))}
# inverse s
S <-diag.inverse(s)
# follow chen's paper to caculate omega
omega<-diag.multi(S, V)
# retrieve omega value from the matrix
omega.diag<-diag(omega )
# summary the omega value
summary(omega.diag)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003435 0.020100 0.038650 0.072460 0.080030 0.920500
# betas.thea<- S %*% Z %*% gamma + (identity-S) %*% betas.hat
# caculate betas.tieda with the formula in Chen's paper
betas.tieda<- omega %*% Z %*% gamma + (identity-omega) %*% betas.hat
# crbetas.tieda<- cromega %*% Z %*% gamma + (identity-cromega) %*% betas.hat
head(betas.tieda)
## [,1]
## [1,] 0.03521062
## [2,] 0.08438259
## [3,] 0.06779637
## [4,] 0.05187560
## [5,] 0.06892069
## [6,] 0.19018058
head(betas.hat)
## [,1]
## [1,] 0.03197222
## [2,] 0.08446667
## [3,] 0.06783732
## [4,] 0.05122222
## [5,] 0.06860287
## [6,] 0.19240670
#regression beta
regbeta <-Z %*% gamma
head(regbeta)
## [,1]
## [1,] 0.08341975
## [2,] 0.07090379
## [3,] 0.06740610
## [4,] 0.08977023
## [5,] 0.08101177
## [6,] 0.12850308
summary(regbeta)
## V1
## Min. :0.06271
## 1st Qu.:0.07491
## Median :0.08667
## Mean :0.10929
## 3rd Qu.:0.11178
## Max. :1.33060
markers1<-as.character(markers)
# combine ensemble_id, betas.hat and betas.tieda
outputVector<-c(markers1,betas.hat,betas.tieda)
write.table(matrix(outputVector,rowLength),file="2016-04-26_hm_tau_hmresults.txt",col.names=FALSE,row.names=FALSE,quote=FALSE)
liver.mouse.eQTL.bayesian<-read.table(file="2016-04-26_hm_tau_hmresults.txt")
colnames(liver.mouse.eQTL.bayesian)<-c( "ensembl_id", "betas.hat","betas.tieda")
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat betas.tieda
## 1 ENSMUSG00000000001 0.03197222 0.03521062
## 2 ENSMUSG00000000049 0.08446667 0.08438259
## 3 ENSMUSG00000000056 0.06783732 0.06779637
## 4 ENSMUSG00000000058 0.05122222 0.05187560
## 5 ENSMUSG00000000085 0.06860287 0.06892069
## 6 ENSMUSG00000000088 0.19240670 0.19018058
# merge dataset with betas.hat and betas.tieda
liver.mouse.eQTL.bayesian.all<- merge(liver.mouse.eQTL.bayesian, merged.mouse.eQTL.min.variance2, by = "ensembl_id")
head(liver.mouse.eQTL.bayesian.all)
## ensembl_id betas.hat betas.tieda lung_pvalue lung.beta
## 1 ENSMUSG00000000001 0.03197222 0.03521062 0.159238340 -0.05980000
## 2 ENSMUSG00000000049 0.08446667 0.08438259 0.896513064 -0.02366071
## 3 ENSMUSG00000000056 0.06783732 0.06779637 0.698124238 -0.01356126
## 4 ENSMUSG00000000058 0.05122222 0.05187560 0.314117648 -0.07813675
## 5 ENSMUSG00000000085 0.06860287 0.06892069 0.363075471 0.05284706
## 6 ENSMUSG00000000088 0.19240670 0.19018058 0.001444641 0.18997619
## lung.beta_se liver_pvalue liver.beta liver.beta_se abs_liver.beta
## 1 0.04174596 5.377718e-01 -0.03197222 0.05124934 0.03197222
## 2 0.18083869 9.084753e-06 0.08446667 0.01561692 0.08446667
## 3 0.03473152 2.986190e-01 0.06783732 0.06405365 0.06783732
## 4 0.07671439 5.846922e-02 -0.05122222 0.02596463 0.05122222
## 5 0.05748721 4.121599e-02 -0.06860287 0.03205858 0.06860287
## 6 0.05579620 1.987922e-05 0.19240670 0.03756628 0.19240670
## abs_lung.beta neg_log_lung_pvalue
## 1 0.05980000 0.79795236
## 2 0.02366071 0.04744338
## 3 0.01356126 0.15606728
## 4 0.07813675 0.50290766
## 5 0.05284706 0.44000309
## 6 0.18997619 2.84024014
write.table(liver.mouse.eQTL.bayesian.all,file="2016-05-09_liver.mouse.eQTL.bayesian.all.txt")
liver.mouse.eQTL.bayesian<-read.table(file="2016-05-09_liver.mouse.eQTL.bayesian.all.txt")
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat betas.tieda lung_pvalue lung.beta
## 1 ENSMUSG00000000001 0.03197222 0.03521062 0.159238340 -0.05980000
## 2 ENSMUSG00000000049 0.08446667 0.08438259 0.896513064 -0.02366071
## 3 ENSMUSG00000000056 0.06783732 0.06779637 0.698124238 -0.01356126
## 4 ENSMUSG00000000058 0.05122222 0.05187560 0.314117648 -0.07813675
## 5 ENSMUSG00000000085 0.06860287 0.06892069 0.363075471 0.05284706
## 6 ENSMUSG00000000088 0.19240670 0.19018058 0.001444641 0.18997619
## lung.beta_se liver_pvalue liver.beta liver.beta_se abs_liver.beta
## 1 0.04174596 5.377718e-01 -0.03197222 0.05124934 0.03197222
## 2 0.18083869 9.084753e-06 0.08446667 0.01561692 0.08446667
## 3 0.03473152 2.986190e-01 0.06783732 0.06405365 0.06783732
## 4 0.07671439 5.846922e-02 -0.05122222 0.02596463 0.05122222
## 5 0.05748721 4.121599e-02 -0.06860287 0.03205858 0.06860287
## 6 0.05579620 1.987922e-05 0.19240670 0.03756628 0.19240670
## abs_lung.beta neg_log_lung_pvalue
## 1 0.05980000 0.79795236
## 2 0.02366071 0.04744338
## 3 0.01356126 0.15606728
## 4 0.07813675 0.50290766
## 5 0.05284706 0.44000309
## 6 0.18997619 2.84024014
liver.mouse.eQTL.bayesian<-subset(liver.mouse.eQTL.bayesian, select = c("ensembl_id", "betas.hat",
"liver.beta_se", "betas.tieda",
"liver_pvalue", "abs_lung.beta",
"abs_lung.beta", "neg_log_lung_pvalue"))
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat liver.beta_se betas.tieda liver_pvalue
## 1 ENSMUSG00000000001 0.03197222 0.05124934 0.03521062 5.377718e-01
## 2 ENSMUSG00000000049 0.08446667 0.01561692 0.08438259 9.084753e-06
## 3 ENSMUSG00000000056 0.06783732 0.06405365 0.06779637 2.986190e-01
## 4 ENSMUSG00000000058 0.05122222 0.02596463 0.05187560 5.846922e-02
## 5 ENSMUSG00000000085 0.06860287 0.03205858 0.06892069 4.121599e-02
## 6 ENSMUSG00000000088 0.19240670 0.03756628 0.19018058 1.987922e-05
## abs_lung.beta abs_lung.beta.1 neg_log_lung_pvalue
## 1 0.05980000 0.05980000 0.79795236
## 2 0.02366071 0.02366071 0.04744338
## 3 0.01356126 0.01356126 0.15606728
## 4 0.07813675 0.07813675 0.50290766
## 5 0.05284706 0.05284706 0.44000309
## 6 0.18997619 0.18997619 2.84024014
# Caculate variance for beta.tieda by following Brian Kulis' lecture notes
# Invert Tau and V
Tau_invert<-diag.inverse(Tau)
V_invert<-diag.inverse(V)
PS_invert<-Tau_invert + V_invert
# PS_invert<-Tau_invert+V_invert%*% Z %*% Z_transpose # previous wrong code
# S in Brian Kulis' lecture note:PS
PS <- diag.inverse(PS_invert)
# retrieve posterior variance
ps<-diag(PS)
range(ps)
## [1] 0.0001343022 0.0359904492
library(reshape)
##
## Attaching package: 'reshape'
## The following objects are masked from 'package:plyr':
##
## rename, round_any
## The following object is masked from 'package:data.table':
##
## melt
# reshape posterior variance to long format
ps.long <- melt(ps)
head(ps.long)
## value
## 1 0.0024611682
## 2 0.0002423764
## 3 0.0037132296
## 4 0.0006627353
## 5 0.0010014297
## 6 0.0013620645
# Caculate sd: square root on variance
ps.long$betas.tieda.se<-(ps.long$value)^0.5
# combine sd to the data.frame
liver.mouse.eQTL.bayesian<-cbind(liver.mouse.eQTL.bayesian,ps.long$betas.tieda.se)
# head(liver.mouse.eQTL.bayesian)
# rename betas.tieda.se
liver.mouse.eQTL.bayesian<-rename(liver.mouse.eQTL.bayesian, c("ps.long$betas.tieda.se"="betas.tieda.se", "liver.beta_se"="betas.hat.se"))
liver.mouse.eQTL.bayesian<-subset(liver.mouse.eQTL.bayesian, select = c("ensembl_id", "betas.hat", "betas.hat.se",
"betas.tieda", "betas.tieda.se",
"liver_pvalue", "abs_lung.beta", "neg_log_lung_pvalue"))
# head(liver.mouse.eQTL.bayesian)
# library(tigerstats)
# pnormGC(0, region="below", mean=0.002352829,sd=0.09972950)
# caculate probability of betas.tieda below 0 based on betas.tieda and standard deviation
liver.mouse.eQTL.bayesian$p.below.0 <- pnorm(0,liver.mouse.eQTL.bayesian$betas.tieda, liver.mouse.eQTL.bayesian$betas.tieda.se)
head(liver.mouse.eQTL.bayesian)
## ensembl_id betas.hat betas.hat.se betas.tieda betas.tieda.se
## 1 ENSMUSG00000000001 0.03197222 0.05124934 0.03521062 0.04961016
## 2 ENSMUSG00000000049 0.08446667 0.01561692 0.08438259 0.01556844
## 3 ENSMUSG00000000056 0.06783732 0.06405365 0.06779637 0.06093627
## 4 ENSMUSG00000000058 0.05122222 0.02596463 0.05187560 0.02574365
## 5 ENSMUSG00000000085 0.06860287 0.03205858 0.06892069 0.03164537
## 6 ENSMUSG00000000088 0.19240670 0.03756628 0.19018058 0.03690616
## liver_pvalue abs_lung.beta neg_log_lung_pvalue p.below.0
## 1 5.377718e-01 0.05980000 0.79795236 2.389308e-01
## 2 9.084753e-06 0.02366071 0.04744338 2.978206e-08
## 3 2.986190e-01 0.01356126 0.15606728 1.329448e-01
## 4 5.846922e-02 0.07813675 0.50290766 2.194794e-02
## 5 4.121599e-02 0.05284706 0.44000309 1.470647e-02
## 6 1.987922e-05 0.18997619 2.84024014 1.281176e-07
dim(liver.mouse.eQTL.bayesian)
## [1] 10601 9
summary(liver.mouse.eQTL.bayesian$betas.tieda.se)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01159 0.02803 0.03887 0.04639 0.05594 0.18970
range(liver.mouse.eQTL.bayesian$p.below.0)
## [1] 0.0000000 0.4855218
write.table(liver.mouse.eQTL.bayesian,file="2016-05-04_liver.mouse.eQTL.bayesian with beta.txt")
liver.mouse.eQTL.bayesian <- read.table(file="2016-05-04_liver.mouse.eQTL.bayesian with beta.txt")
# head(liver.mouse.eQTL.bayesian)
# summary(liver.mouse.eQTL.bayesian$liver_residual_variance)
liver.mouse.eQTL.bayesian.4tau <- liver.mouse.eQTL.bayesian
# colnames(liver.mouse.eQTL.bayesian.4tau) <- c("ensembl_id", "betas.hat", "betas.hat.se", "betas.tieda", "betas.tieda.se", "liver_residual_variance", "liver_pvalue", "abs_lung.beta",
# "neg_log_lung_pvalue", "p.below.0", "betas.tieda2m", "betas.tieda3rd", "betas.tieda4max")
# head(liver.mouse.eQTL.bayesian.4tau)
# Introduce weight (Tmm) to adjust Tau with neg_log_lung_pvalue
liver.mouse.eQTL.bayesian.4tau$fzm <- liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue
# caculate ratio_fzm
liver.mouse.eQTL.bayesian.4tau$ratio_fzm <- max(liver.mouse.eQTL.bayesian.4tau$fzm)/liver.mouse.eQTL.bayesian.4tau$fzm
range(liver.mouse.eQTL.bayesian.4tau$ratio_fzm)
## [1] 1 Inf
# set up threshold for ratio_fzm and caculate updated ratio_fzm (nratio_fzm)
threshold <- 0.05
liver.mouse.eQTL.bayesian.4tau$nratio_fzm <- liver.mouse.eQTL.bayesian.4tau$ratio_fzm
liver.mouse.eQTL.bayesian.4tau$nratio_fzm[liver.mouse.eQTL.bayesian.4tau$ratio_fzm > max(liver.mouse.eQTL.bayesian.4tau$fzm)/(-log10(threshold))] <- max(liver.mouse.eQTL.bayesian.4tau$fzm)/(-log10(threshold))
# liver.mouse.eQTL.bayesian.4tau$nratio_fzm[liver.mouse.eQTL.bayesian.4tau$nratio_fzm >= max(liver.mouse.eQTL.bayesian.4tau$fzm)/(-log(threshold))] <- liver.mouse.eQTL.bayesian.4tau$ratio_fzm
summary(liver.mouse.eQTL.bayesian.4tau$nratio_fzm)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 25.04 35.84 29.41 35.84 35.84
# compare bayesian prediction to the traditional method
# evaluate the predition with alle specific expreesion in the liver: Sandrine Lagarrigue's paper
liver.ASE <- read.csv(file= "ASE.genetics.113.153882-6.csv")
dim(liver.ASE)
## [1] 1191 19
head(liver.ASE)
## replicate chr startExon geneID SNPperExon3
## 1 M.CH. DxB and BxD 1 9535488 Rrs1 15
## 2 M.HF DxB and BxD 1 9535488 Rrs1 14
## 3 F.HF DxB and BxD 1 9535488 Rrs1 15
## 4 M.CH. DxB and BxD 1 37473929 6330578E17Rik 2
## 5 M.CH. DxB and BxD 1 58169979 Aox3 3
## 6 M.HF DxB and BxD 1 58169979 Aox3 3
## sumBperExon.DxB4 sumDperExon.DxB4 sumBperExon.BxD4 sumDperExon.BxD4
## 1 45 19 50 25
## 2 74 39 66 30
## 3 76 20 77 40
## 4 78 32 47 17
## 5 473 82 225 27
## 6 252 56 263 53
## FCadd1.DxB5 FCadd1.BxD5 BonBD.DxB6 BonBD.BxD6 pvalBH.DxB7 pvalBH.BxD7
## 1 2.30 1.96 0.70 0.67 1.0e-02 3.7e-02
## 2 1.88 2.16 0.65 0.69 1.1e-02 3.3e-03
## 3 3.67 1.90 0.79 0.66 4.9e-07 5.1e-03
## 4 2.39 2.67 0.71 0.73 1.9e-04 2.2e-03
## 5 5.71 8.07 0.85 0.89 3.6e-64 5.4e-37
## 6 4.44 4.89 0.82 0.83 8.7e-28 1.1e-31
## UTR5 UTR3 strand exonCount
## 1 0 0 + 1
## 2 0 0 + 1
## 3 0 0 + 1
## 4 0 0 - 3
## 5 0 0 + 35
## 6 0 0 + 35
# 440 unique gene ID
length(unique(liver.ASE$geneID))
## [1] 440
# verify ASE table
liver.ASE1 <- liver.ASE[which(liver.ASE$replicate == "M.CH. DxB and BxD"), ]
liver.ASE2 <- liver.ASE[which(liver.ASE$replicate == "M.HF DxB and BxD"), ]
liver.ASE3 <- liver.ASE[which(liver.ASE$replicate == "F.HF DxB and BxD"), ]
length(unique(liver.ASE1$geneID))
## [1] 272
length(unique(liver.ASE2$geneID))
## [1] 275
length(unique(liver.ASE3$geneID))
## [1] 304
(length(unique(liver.ASE1$geneID))+length(unique(liver.ASE2$geneID))+length(unique(liver.ASE3$geneID)))/3
## [1] 283.6667
# As claimed in the paper: averaged 284 ASE for each replicate
sub.liver.ASE <-liver.ASE1
summary(sub.liver.ASE$pvalBH.DxB7)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000000 0.0000000 0.0000058 0.0084070 0.0031000 0.1000000
sub.liver.ASE1 <- subset(sub.liver.ASE, pvalBH.DxB7 < 0.00000000000001)
sub.liver.ASE2 <- subset(sub.liver.ASE, pvalBH.DxB7 >= 0.00000000000001 & pvalBH.DxB7 < 0.0000058)
sub.liver.ASE3 <- subset(sub.liver.ASE, pvalBH.DxB7 >= 0.0000058 & pvalBH.DxB7 < 0.0031000)
sub.liver.ASE4 <- subset(sub.liver.ASE, pvalBH.DxB7 >= 0.0031000)
dim(sub.liver.ASE1)
## [1] 89 19
dim(sub.liver.ASE2)
## [1] 97 19
dim(sub.liver.ASE3)
## [1] 93 19
dim(sub.liver.ASE4)
## [1] 94 19
# Subset liver ASE with different conditions
# sub.liver.ASE <-liver.ASE[which(liver.ASE$pvalBH.DxB7 < 0.05 & liver.ASE$pvalBH.BxD7 < 0.05), ]
# sub.liver.ASE <- sub.liver.ASE[order(sub.liver.ASE$pvalBH.BxD7), ]
# summary(sub.liver.ASE$pvalBH.BxD7)
# sub.liver.ASE <- sub.liver.ASE[which(sub.liver.ASE$pvalBH.DxB7 <= 9.0e-10 ), ]
# sub.liver.ASE <- sub.liver.ASE[which(sub.liver.ASE$pvalBH.BxD7 <= 9.0e-10 ), ]
# sub.liver.ASE <- sub.liver.ASE[ sub.liver.ASE$geneID %in% names(table(sub.liver.ASE$geneID))[table(sub.liver.ASE$geneID) >1] , ]
# check the remain gene number after subsetting
dim(sub.liver.ASE)
## [1] 373 19
liver.ASE.symbol <- unique(sub.liver.ASE$geneID)
liver.ASE.symbol1 <- unique(sub.liver.ASE1$geneID)
liver.ASE.symbol2 <- unique(sub.liver.ASE2$geneID)
liver.ASE.symbol3 <- unique(sub.liver.ASE3$geneID)
liver.ASE.symbol4 <- unique(sub.liver.ASE4$geneID)
length(liver.ASE.symbol)
## [1] 272
# Annoate gene symbol wiht ensemble.ID
library(biomaRt)
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
liver.ASE.ensembl <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol", values =liver.ASE.symbol, mart=mouse)
liver.ASE.ensembl1 <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol", values =liver.ASE.symbol1, mart=mouse)
liver.ASE.ensembl2 <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol", values =liver.ASE.symbol2, mart=mouse)
liver.ASE.ensembl3 <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol", values =liver.ASE.symbol3, mart=mouse)
liver.ASE.ensembl4 <- getBM( attributes=c("ensembl_gene_id", "mgi_symbol") , filters=
"mgi_symbol", values =liver.ASE.symbol4, mart=mouse)
dim(liver.ASE.ensembl)
## [1] 241 2
liver.ASE.ensembl <- unique(liver.ASE.ensembl)
# delete liver ASE ensemble ID which are not in the liver.mouse.eQTL.bayesian data frame
liver.ASE.ensembl <- liver.ASE.ensembl[liver.ASE.ensembl$ensembl_gene_id %in% liver.mouse.eQTL.bayesian.4tau$ensembl_id, ]
dim(liver.ASE.ensembl)
## [1] 195 2
# create indicator for ASE true or not
# liver.mouse.eQTL.bayesian.4tau$eqtl[liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl1$ensembl_gene_id] <- 1
# liver.mouse.eQTL.bayesian.4tau$eqtl[liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl2$ensembl_gene_id] <- 2
# liver.mouse.eQTL.bayesian.4tau$eqtl[liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl3$ensembl_gene_id] <- 3
# liver.mouse.eQTL.bayesian.4tau$eqtl[liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl4$ensembl_gene_id] <- 4
# liver.mouse.eQTL.bayesian.4tau$eqtl[!liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id] <- 5
liver.mouse.eQTL.bayesian.4tau$eqtl[liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id] <- 1
liver.mouse.eQTL.bayesian.4tau$eqtl[!liver.mouse.eQTL.bayesian.4tau$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id] <- 0
summary(liver.mouse.eQTL.bayesian.4tau$eqtl)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01839 0.00000 1.00000
liver.mouse.eQTL.bayesian.4tau$neg_log_liver_pvalue <- -log10(liver.mouse.eQTL.bayesian.4tau$liver_pvalue)
head(liver.mouse.eQTL.bayesian.4tau)
## ensembl_id betas.hat betas.hat.se betas.tieda betas.tieda.se
## 1 ENSMUSG00000000001 0.03197222 0.05124934 0.03521062 0.04961016
## 2 ENSMUSG00000000049 0.08446667 0.01561692 0.08438259 0.01556844
## 3 ENSMUSG00000000056 0.06783732 0.06405365 0.06779637 0.06093627
## 4 ENSMUSG00000000058 0.05122222 0.02596463 0.05187560 0.02574365
## 5 ENSMUSG00000000085 0.06860287 0.03205858 0.06892069 0.03164537
## 6 ENSMUSG00000000088 0.19240670 0.03756628 0.19018058 0.03690616
## liver_pvalue abs_lung.beta neg_log_lung_pvalue p.below.0 fzm
## 1 5.377718e-01 0.05980000 0.79795236 2.389308e-01 0.79795236
## 2 9.084753e-06 0.02366071 0.04744338 2.978206e-08 0.04744338
## 3 2.986190e-01 0.01356126 0.15606728 1.329448e-01 0.15606728
## 4 5.846922e-02 0.07813675 0.50290766 2.194794e-02 0.50290766
## 5 4.121599e-02 0.05284706 0.44000309 1.470647e-02 0.44000309
## 6 1.987922e-05 0.18997619 2.84024014 1.281176e-07 2.84024014
## ratio_fzm nratio_fzm eqtl neg_log_liver_pvalue
## 1 58.43301 35.83834 0 0.2694020
## 2 982.78748 35.83834 0 5.0416869
## 3 298.76061 35.83834 0 0.5248826
## 4 92.71435 35.83834 0 1.2330727
## 5 105.96916 35.83834 0 1.3849342
## 6 16.41648 16.41648 0 4.7016006
by(liver.mouse.eQTL.bayesian.4tau[, c(1, 7, 9, 14)], liver.mouse.eQTL.bayesian.4tau[, "eqtl"], summary)
## liver.mouse.eQTL.bayesian.4tau[, "eqtl"]: 0
## ensembl_id abs_lung.beta p.below.0
## ENSMUSG00000000001: 1 Min. :0.00000 Min. :0.00000
## ENSMUSG00000000049: 1 1st Qu.:0.03487 1st Qu.:0.01514
## ENSMUSG00000000056: 1 Median :0.06815 Median :0.09871
## ENSMUSG00000000058: 1 Mean :0.13143 Mean :0.13682
## ENSMUSG00000000085: 1 3rd Qu.:0.13975 3rd Qu.:0.23055
## ENSMUSG00000000088: 1 Max. :3.66099 Max. :0.48552
## (Other) :10400
## neg_log_liver_pvalue
## Min. : 0.0000
## 1st Qu.: 0.2853
## Median : 0.6341
## Mean : 1.3529
## 3rd Qu.: 1.3629
## Max. :33.3106
##
## --------------------------------------------------------
## liver.mouse.eQTL.bayesian.4tau[, "eqtl"]: 1
## ensembl_id abs_lung.beta p.below.0
## ENSMUSG00000000275: 1 Min. :0.001944 Min. :0.0000000
## ENSMUSG00000000673: 1 1st Qu.:0.076576 1st Qu.:0.0000000
## ENSMUSG00000001467: 1 Median :0.164283 Median :0.0000028
## ENSMUSG00000001473: 1 Mean :0.299202 Mean :0.0376244
## ENSMUSG00000001604: 1 3rd Qu.:0.360720 3rd Qu.:0.0194118
## ENSMUSG00000002395: 1 Max. :3.659845 Max. :0.4351068
## (Other) :189
## neg_log_liver_pvalue
## Min. : 0.00735
## 1st Qu.: 1.29820
## Median : 4.36302
## Mean : 5.73612
## 3rd Qu.: 9.29962
## Max. :23.22672
##
library(ggplot2)
boxplot(neg_log_liver_pvalue ~ eqtl,data=liver.mouse.eQTL.bayesian.4tau, main="liver.mouse.eQTL",
xlab="ASE cutoff by p value", ylab="liver neg log p")
boxplot(neg_log_lung_pvalue ~ eqtl,data=liver.mouse.eQTL.bayesian.4tau, main="lung.mouse.eQTL",
xlab="ASE cutoff by p value", ylab="lung neg log p")
liver.mouse.eQTL.bayesian.4tau.ase <- liver.mouse.eQTL.bayesian.4tau[liver.mouse.eQTL.bayesian.4tau$eqtl == 1, ]
plot(liver.mouse.eQTL.bayesian.4tau$neg_log_liver_pvalue, liver.mouse.eQTL.bayesian.4tau$neg_log_lung_pvalue, col=factor(liver.mouse.eQTL.bayesian.4tau$eqtl), xlab="neg_log_liver_pvalue", ylab="neg_log_lung_pvalue" )
legend("topright", cex = .75, inset=.05, c("ASE","others"), text.col = c("red", "black"), horiz=TRUE)
plot(liver.mouse.eQTL.bayesian.4tau.ase$neg_log_liver_pvalue, liver.mouse.eQTL.bayesian.4tau.ase$neg_log_lung_pvalue, , col="red", xlab="neg_log_liver_pvalue", ylab="neg_log_lung_pvalue")
legend("topright", cex = .75, inset=.05, c("ASE"), text.col = c("red"), horiz=TRUE)
plot(liver.mouse.eQTL.bayesian.4tau$betas.hat, liver.mouse.eQTL.bayesian.4tau$betas.tieda, col=factor(liver.mouse.eQTL.bayesian.4tau$eqtl) )
legend("topright", cex = .75, inset=.05, c("ASE","others"), text.col = c("red", "black"), horiz=TRUE)
plot(liver.mouse.eQTL.bayesian.4tau.ase$betas.hat, liver.mouse.eQTL.bayesian.4tau.ase$betas.tieda, col="red")
legend("topright", cex = .75, inset=.05, c("ASE"), text.col = c("red"), horiz=TRUE)
cor(liver.mouse.eQTL.bayesian.4tau.ase$neg_log_liver_pvalue, liver.mouse.eQTL.bayesian.4tau.ase$neg_log_lung_pvalue)
## [1] 0.4483848
length(liver.mouse.eQTL.bayesian.4tau.ase$neg_log_liver_pvalue)
## [1] 195
plot(liver.mouse.eQTL.bayesian.4tau$neg_log_liver_pvalue, liver.mouse.eQTL.bayesian.4tau$p.below.0, col=factor(liver.mouse.eQTL.bayesian.4tau$eqtl) )
legend("topright", cex = .75, inset=.05, c("ASE","others"), text.col = c("red", "black"), horiz=TRUE)
plot(liver.mouse.eQTL.bayesian.4tau.ase$neg_log_liver_pvalue, liver.mouse.eQTL.bayesian.4tau.ase$p.below.0, col="red")
legend("topright", cex = .75, inset=.05, c("ASE"), text.col = c("red"), horiz=TRUE)
head(liver.mouse.eQTL.bayesian.4tau)
## ensembl_id betas.hat betas.hat.se betas.tieda betas.tieda.se
## 1 ENSMUSG00000000001 0.03197222 0.05124934 0.03521062 0.04961016
## 2 ENSMUSG00000000049 0.08446667 0.01561692 0.08438259 0.01556844
## 3 ENSMUSG00000000056 0.06783732 0.06405365 0.06779637 0.06093627
## 4 ENSMUSG00000000058 0.05122222 0.02596463 0.05187560 0.02574365
## 5 ENSMUSG00000000085 0.06860287 0.03205858 0.06892069 0.03164537
## 6 ENSMUSG00000000088 0.19240670 0.03756628 0.19018058 0.03690616
## liver_pvalue abs_lung.beta neg_log_lung_pvalue p.below.0 fzm
## 1 5.377718e-01 0.05980000 0.79795236 2.389308e-01 0.79795236
## 2 9.084753e-06 0.02366071 0.04744338 2.978206e-08 0.04744338
## 3 2.986190e-01 0.01356126 0.15606728 1.329448e-01 0.15606728
## 4 5.846922e-02 0.07813675 0.50290766 2.194794e-02 0.50290766
## 5 4.121599e-02 0.05284706 0.44000309 1.470647e-02 0.44000309
## 6 1.987922e-05 0.18997619 2.84024014 1.281176e-07 2.84024014
## ratio_fzm nratio_fzm eqtl neg_log_liver_pvalue
## 1 58.43301 35.83834 0 0.2694020
## 2 982.78748 35.83834 0 5.0416869
## 3 298.76061 35.83834 0 0.5248826
## 4 92.71435 35.83834 0 1.2330727
## 5 105.96916 35.83834 0 1.3849342
## 6 16.41648 16.41648 0 4.7016006
# Optimizing rho and adjust the weight
library(reshape)
rho.optimization <- matrix(0, nrow=nrow(liver.mouse.eQTL.bayesian.4tau), ncol=8)
colnames(rho.optimization)<-c("class", "rho","tmm","tau", "omega","beta_tieda", "n.betas.tieda.se","p.below.0")
nomega.diag<-diag(omega )
rhoclass <- seq(1,11, by=2)
rho <- rhoclass*tau
result <- NULL
for (i in 1:length(rho)) {
rho.optimization[ ,1] <- rhoclass[i]
rho.optimization[ ,2] <- rho[i]
rho.optimization[ ,3] <- (rho[i]/tau)^liver.mouse.eQTL.bayesian.4tau$nratio_fzm
rho.optimization[ ,4] <-tau*((rho[i]/tau)^liver.mouse.eQTL.bayesian.4tau$nratio_fzm)
nTau<- diag(rho.optimization[ ,3], rowLength, rowLength)
ns<-V + nTau
nS <- diag.inverse(ns)
nomega<-diag.multi(nS, V)
# nomega <- diag(0, rowLength, rowLength) # set nomega to 0 for code checking
# nomega <- diag(1, rowLength, rowLength) # set nomega to 1 for code checking
rho.optimization[ ,5] <- diag(nomega )
rho.optimization[ ,6] <- nomega %*% Z %*% gamma + (identity-nomega) %*% betas.hat
nTau_invert<-diag.inverse(nTau)
V_invert<-diag.inverse(V)
nPS_invert<-nTau_invert+ V_invert
# nPS_invert<-nTau_invert+ diag.multi(diag.multi(V_invert, Z_transpose), Z) # previous wrong code
nPS<-diag.inverse(nPS_invert)
nps<-diag(nPS)
nps.long <- melt(nps)
rho.optimization[ ,7] <-(nps.long$value)^0.5
rho.optimization[ ,8] <- pnorm(0, rho.optimization[ ,6], rho.optimization[ ,7])
result <- rbind(result,rho.optimization)
}
dim(result)
## [1] 63606 8
head(result)
## class rho tmm tau omega beta_tieda
## [1,] 1 0.03909984 1 0.03909984 0.0026196144 0.03210699
## [2,] 1 0.03909984 1 0.03909984 0.0002438288 0.08446336
## [3,] 1 0.03909984 1 0.03909984 0.0040861058 0.06783556
## [4,] 1 0.03909984 1 0.03909984 0.0006737081 0.05124819
## [5,] 1 0.03909984 1 0.03909984 0.0010266974 0.06861561
## [6,] 1 0.03909984 1 0.03909984 0.0014092365 0.19231664
## n.betas.tieda.se p.below.0
## [1,] 0.05118217 2.652286e-01
## [2,] 0.01561502 3.166925e-08
## [3,] 0.06392265 1.442965e-01
## [4,] 0.02595589 2.416619e-02
## [5,] 0.03204212 1.612012e-02
## [6,] 0.03753980 1.503509e-07
write.table(result, file="2016-05-04_liver.mouse.eQTL.bayesian.result.txt",col.names=TRUE,row.names=FALSE,quote=FALSE)
liver.mouse.eQTL.bayesian.result <- read.table(file="2016-05-04_liver.mouse.eQTL.bayesian.result.txt", header=T)
result.df <-liver.mouse.eQTL.bayesian.result
result.df$class <- factor(result.df$class)
# combine liver.mouse.eqtl.bayesian and rho.optimization.result for ploting
a <-liver.mouse.eQTL.bayesian.4tau[, c(1:2, 6, 7)]
a <-rbind(a, a, a, a, a, a)
dim(a)
## [1] 63606 4
new.result.df<-cbind(a, result.df)
head(new.result.df)
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 1 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 1
## 2 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 1
## 3 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 1
## 4 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 1
## 5 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 1
## 6 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 1
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 1 0.03909984 1 0.03909984 0.0026196144 0.03210699 0.05118217
## 2 0.03909984 1 0.03909984 0.0002438288 0.08446336 0.01561502
## 3 0.03909984 1 0.03909984 0.0040861058 0.06783556 0.06392265
## 4 0.03909984 1 0.03909984 0.0006737081 0.05124819 0.02595589
## 5 0.03909984 1 0.03909984 0.0010266974 0.06861561 0.03204212
## 6 0.03909984 1 0.03909984 0.0014092365 0.19231664 0.03753980
## p.below.0
## 1 2.652286e-01
## 2 3.166925e-08
## 3 1.442965e-01
## 4 2.416619e-02
## 5 1.612012e-02
## 6 1.503509e-07
new.result.df2 <- new.result.df
by(new.result.df2, new.result.df2[, "class"], head)
## new.result.df2[, "class"]: 1
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 1 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 1
## 2 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 1
## 3 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 1
## 4 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 1
## 5 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 1
## 6 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 1
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 1 0.03909984 1 0.03909984 0.0026196144 0.03210699 0.05118217
## 2 0.03909984 1 0.03909984 0.0002438288 0.08446336 0.01561502
## 3 0.03909984 1 0.03909984 0.0040861058 0.06783556 0.06392265
## 4 0.03909984 1 0.03909984 0.0006737081 0.05124819 0.02595589
## 5 0.03909984 1 0.03909984 0.0010266974 0.06861561 0.03204212
## 6 0.03909984 1 0.03909984 0.0014092365 0.19231664 0.03753980
## p.below.0
## 1 2.652286e-01
## 2 3.166925e-08
## 3 1.442965e-01
## 4 2.416619e-02
## 5 1.612012e-02
## 6 1.503509e-07
## --------------------------------------------------------
## new.result.df2[, "class"]: 3
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 11000 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 3
## 21000 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 3
## 31000 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 3
## 41000 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 3
## 51000 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 3
## 61000 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 3
## rho tmm tau omega beta_tieda
## 11000 0.1172995 1.256708e+17 4.913710e+15 2.089980e-20 0.03197222
## 21000 0.1172995 1.256708e+17 4.913710e+15 1.940691e-21 0.08446667
## 31000 0.1172995 1.256708e+17 4.913710e+15 3.264775e-20 0.06783732
## 41000 0.1172995 1.256708e+17 4.913710e+15 5.364508e-21 0.05122222
## 51000 0.1172995 1.256708e+17 4.913710e+15 8.178131e-21 0.06860287
## 61000 0.1172995 6.802265e+07 2.659675e+06 2.074640e-11 0.19240670
## n.betas.tieda.se p.below.0
## 11000 0.05124934 2.663610e-01
## 21000 0.01561692 3.174851e-08
## 31000 0.06405365 1.447839e-01
## 41000 0.02596463 2.426095e-02
## 51000 0.03205858 1.618052e-02
## 61000 0.03756628 1.513221e-07
## --------------------------------------------------------
## new.result.df2[, "class"]: 5
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 11002 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 5
## 21002 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 5
## 31002 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 5
## 41002 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 5
## 51002 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 5
## 61002 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 5
## rho tmm tau omega beta_tieda
## 11002 0.1954992 1.121827e+25 4.386325e+23 2.341266e-28 0.03197222
## 21002 0.1954992 1.121827e+25 4.386325e+23 2.174028e-29 0.08446667
## 31002 0.1954992 1.121827e+25 4.386325e+23 3.657312e-28 0.06783732
## 41002 0.1954992 1.121827e+25 4.386325e+23 6.009504e-29 0.05122222
## 51002 0.1954992 1.121827e+25 4.386325e+23 9.161420e-29 0.06860287
## 61002 0.1954992 2.982841e+11 1.166286e+10 4.731146e-15 0.19240670
## n.betas.tieda.se p.below.0
## 11002 0.05124934 2.663610e-01
## 21002 0.01561692 3.174851e-08
## 31002 0.06405365 1.447839e-01
## 41002 0.02596463 2.426095e-02
## 51002 0.03205858 1.618052e-02
## 61002 0.03756628 1.513221e-07
## --------------------------------------------------------
## new.result.df2[, "class"]: 7
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 11004 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 7
## 21004 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 7
## 31004 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 7
## 41004 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 7
## 51004 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 7
## 61004 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 7
## rho tmm tau omega beta_tieda
## 11004 0.2736989 1.936031e+30 7.569850e+28 1.356639e-33 0.03197222
## 21004 0.2736989 1.936031e+30 7.569850e+28 1.259733e-34 0.08446667
## 31004 0.2736989 1.936031e+30 7.569850e+28 2.119218e-33 0.06783732
## 41004 0.2736989 1.936031e+30 7.569850e+28 3.482188e-34 0.05122222
## 51004 0.2736989 1.936031e+30 7.569850e+28 5.308555e-34 0.06860287
## 61004 0.2736989 7.473751e+13 2.922225e+12 1.888242e-17 0.19240670
## n.betas.tieda.se p.below.0
## 11004 0.05124934 2.663610e-01
## 21004 0.01561692 3.174851e-08
## 31004 0.06405365 1.447839e-01
## 41004 0.02596463 2.426095e-02
## 51004 0.03205858 1.618052e-02
## 61004 0.03756628 1.513221e-07
## --------------------------------------------------------
## new.result.df2[, "class"]: 9
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 11006 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 9
## 21006 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 9
## 31006 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 9
## 41006 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 9
## 51006 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 9
## 61006 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 9
## rho tmm tau omega beta_tieda
## 11006 0.3518986 1.579316e+34 6.175101e+32 1.663059e-37 0.03197222
## 21006 0.3518986 1.579316e+34 6.175101e+32 1.544265e-38 0.08446667
## 31006 0.3518986 1.579316e+34 6.175101e+32 2.597878e-37 0.06783732
## 41006 0.3518986 1.579316e+34 6.175101e+32 4.268698e-38 0.05122222
## 51006 0.3518986 1.579316e+34 6.175101e+32 6.507580e-38 0.06860287
## 61006 0.3518986 4.627081e+15 1.809181e+14 3.049926e-19 0.19240670
## n.betas.tieda.se p.below.0
## 11006 0.05124934 2.663610e-01
## 21006 0.01561692 3.174851e-08
## 31006 0.06405365 1.447839e-01
## 41006 0.02596463 2.426095e-02
## 51006 0.03205858 1.618052e-02
## 61006 0.03756628 1.513221e-07
## --------------------------------------------------------
## new.result.df2[, "class"]: 11
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 11008 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 11
## 21008 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 11
## 31008 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 11
## 41008 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 11
## 51008 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 11
## 61008 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 11
## rho tmm tau omega beta_tieda
## 11008 0.4300983 2.097910e+37 8.202794e+35 1.251958e-40 0.03197222
## 21008 0.4300983 2.097910e+37 8.202794e+35 1.162530e-41 0.08446667
## 31008 0.4300983 2.097910e+37 8.202794e+35 1.955695e-40 0.06783732
## 41008 0.4300983 2.097910e+37 8.202794e+35 3.213495e-41 0.05122222
## 51008 0.4300983 2.097910e+37 8.202794e+35 4.898936e-41 0.06860287
## 61008 0.4300983 1.247402e+17 4.877324e+15 1.131331e-20 0.19240670
## n.betas.tieda.se p.below.0
## 11008 0.05124934 2.663610e-01
## 21008 0.01561692 3.174851e-08
## 31008 0.06405365 1.447839e-01
## 41008 0.02596463 2.426095e-02
## 51008 0.03205858 1.618052e-02
## 61008 0.03756628 1.513221e-07
# choose different rho class for plotting
new.result.df2.rho1 <- new.result.df2[new.result.df2$class == 1, ]
head(new.result.df2.rho1)
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 1 ENSMUSG00000000001 0.03197222 5.377718e-01 0.05980000 1
## 2 ENSMUSG00000000049 0.08446667 9.084753e-06 0.02366071 1
## 3 ENSMUSG00000000056 0.06783732 2.986190e-01 0.01356126 1
## 4 ENSMUSG00000000058 0.05122222 5.846922e-02 0.07813675 1
## 5 ENSMUSG00000000085 0.06860287 4.121599e-02 0.05284706 1
## 6 ENSMUSG00000000088 0.19240670 1.987922e-05 0.18997619 1
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 1 0.03909984 1 0.03909984 0.0026196144 0.03210699 0.05118217
## 2 0.03909984 1 0.03909984 0.0002438288 0.08446336 0.01561502
## 3 0.03909984 1 0.03909984 0.0040861058 0.06783556 0.06392265
## 4 0.03909984 1 0.03909984 0.0006737081 0.05124819 0.02595589
## 5 0.03909984 1 0.03909984 0.0010266974 0.06861561 0.03204212
## 6 0.03909984 1 0.03909984 0.0014092365 0.19231664 0.03753980
## p.below.0
## 1 2.652286e-01
## 2 3.166925e-08
## 3 1.442965e-01
## 4 2.416619e-02
## 5 1.612012e-02
## 6 1.503509e-07
tail(new.result.df2.rho1)
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 10596 ENSMUSG00000099041 0.12408929 0.02861902 0.10059091 1
## 10597 ENSMUSG00000099083 0.06022624 0.05474021 0.12891000 1
## 10598 ENSMUSG00000099116 0.01886111 0.50556296 0.03457738 1
## 10599 ENSMUSG00000099164 0.07874163 0.10653819 0.08603846 1
## 10600 ENSMUSG00000099262 0.01569378 0.69418700 0.18600427 1
## 10601 ENSMUSG00000099305 0.08752778 0.13257346 0.06246296 1
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 10596 0.03909984 1 0.03909984 0.0028830128 0.12401276 0.05369369
## 10597 0.03909984 1 0.03909984 0.0009016505 0.06026874 0.03002750
## 10598 0.03909984 1 0.03909984 0.0007814649 0.01890474 0.02795469
## 10599 0.03909984 1 0.03909984 0.0022245669 0.07877225 0.04716531
## 10600 0.03909984 1 0.03909984 0.0015582080 0.01586742 0.03947414
## 10601 0.03909984 1 0.03909984 0.0031821105 0.08751764 0.05641020
## p.below.0
## 10596 0.01045422
## 10597 0.02236853
## 10598 0.24943672
## 10599 0.04744674
## 10600 0.34385310
## 10601 0.06039686
# rank with p.below.0: bayesian modeling
new.result.df2.rho1$rank.p.below.0 <- rank(new.result.df2.rho1$p.below.0)
# rank with liver p value: traditionsl linear regression (one step regression)
new.result.df2.rho1$rank.liver.pvalue <- rank(new.result.df2.rho1$liver_pvalue)
new.result.df2.rho1 <- new.result.df2.rho1[order(new.result.df2.rho1$rank.p.below.0), ]
head(new.result.df2.rho1)
## ensembl_id betas.hat liver_pvalue abs_lung.beta class
## 543 ENSMUSG00000004552 5.183400 1.338603e-25 2.0537260 1
## 2577 ENSMUSG00000022680 4.190552 4.891311e-34 2.7728167 1
## 9323 ENSMUSG00000057132 2.458249 7.377910e-27 3.0932000 1
## 612 ENSMUSG00000005161 1.978946 5.933016e-24 1.9872891 1
## 5735 ENSMUSG00000031601 1.635389 2.287024e-22 0.8202054 1
## 5714 ENSMUSG00000031556 1.475527 8.103784e-22 1.1534203 1
## rho tmm tau omega beta_tieda n.betas.tieda.se
## 543 0.03909984 1 0.03909984 0.018383399 5.102340 0.13558539
## 2577 0.03909984 1 0.03909984 0.003004199 4.181036 0.05481057
## 9323 0.03909984 1 0.03909984 0.003400866 2.453745 0.05831695
## 612 0.03909984 1 0.03909984 0.003587004 1.974541 0.05989160
## 5735 0.03909984 1 0.03909984 0.003204701 1.631259 0.05661008
## 5714 0.03909984 1 0.03909984 0.002865289 1.472623 0.05352840
## p.below.0 rank.p.below.0 rank.liver.pvalue
## 543 0.000000e+00 2 3
## 2577 0.000000e+00 2 1
## 9323 0.000000e+00 2 2
## 612 1.145955e-238 4 4
## 5735 6.818177e-183 5 5
## 5714 6.472052e-167 6 6
# caculate TPR: true positive rate
# caculate PPV: positive predictive rate
result.rho1 <- matrix(, nrow(new.result.df2.rho1), 8)
colnames(result.rho1)<-c("bayrank","bayppv","bay_TPR","bay_FPR", "orirank","orippv","ori_TPR","ori_FPR" )
for (i in 1:nrow(new.result.df2.rho1))
{
newdata1.rho1 <- subset(new.result.df2.rho1, rank.p.below.0 <= i)
overlap.newdata1.rho1 <- newdata1.rho1[newdata1.rho1$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
result.rho1[i, 1] <- i
result.rho1[i, 2] <- nrow(overlap.newdata1.rho1)/nrow(newdata1.rho1)
result.rho1[i, 3] <- nrow(overlap.newdata1.rho1)/nrow(liver.ASE.ensembl)
newdata2.rho1 <- subset(new.result.df2.rho1, rank.liver.pvalue <= i)
overlap.newdata2.rho1 <- newdata2.rho1[newdata2.rho1$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
result.rho1[i, 5] <- i
result.rho1[i, 6] <- nrow(overlap.newdata2.rho1)/nrow(newdata2.rho1)
result.rho1[i, 7] <- nrow(overlap.newdata2.rho1)/nrow(liver.ASE.ensembl)
}
head(result.rho1)
## bayrank bayppv bay_TPR bay_FPR orirank orippv ori_TPR
## [1,] 1 NaN 0.000000000 NA 1 0.0000000 0.000000000
## [2,] 2 0.0000000 0.000000000 NA 2 0.0000000 0.000000000
## [3,] 3 0.0000000 0.000000000 NA 3 0.0000000 0.000000000
## [4,] 4 0.2500000 0.005128205 NA 4 0.2500000 0.005128205
## [5,] 5 0.2000000 0.005128205 NA 5 0.2000000 0.005128205
## [6,] 6 0.1666667 0.005128205 NA 6 0.1666667 0.005128205
## ori_FPR
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NA
tail(result.rho1)
## bayrank bayppv bay_TPR bay_FPR orirank orippv ori_TPR
## [10596,] 10596 0.01840317 1 NA 10596 0.01840317 1
## [10597,] 10597 0.01840143 1 NA 10597 0.01840143 1
## [10598,] 10598 0.01839970 1 NA 10598 0.01839970 1
## [10599,] 10599 0.01839796 1 NA 10599 0.01839796 1
## [10600,] 10600 0.01839623 1 NA 10600 0.01839623 1
## [10601,] 10601 0.01839449 1 NA 10601 0.01839449 1
## ori_FPR
## [10596,] NA
## [10597,] NA
## [10598,] NA
## [10599,] NA
## [10600,] NA
## [10601,] NA
# ploting "True positive rate"
plot(result.rho1[, 1], result.rho1[, 3], type="l", col="red", xlab="Ranking", ylab="TPR", ylim=c(0, 1) )
par(new=TRUE)
plot( result.rho1[, 1], result.rho1[, 7], type="l", col="green", xlab="Ranking", ylab="TPR", ylim=c(0, 1))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
plot(result.rho1[, 1], result.rho1[, 3], type="l", col="red", xlab="Ranking", ylab="TPR", ylim=c(0, 0.4) , xlim=c(0, 300))
par(new=TRUE)
plot( result.rho1[, 1], result.rho1[, 7], type="l", col="green", xlab="Ranking", ylab="TPR", ylim=c(0, 0.4), xlim=c(0, 300))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
# ploting "positive predictive value"
plot(result.rho1[, 1], result.rho1[, 2], type="l", col="red", xlab="Ranking", ylab="PPV", ylim=c(0, 1))
par(new=TRUE)
plot(result.rho1[, 5], result.rho1[, 6], type="l", col="green", xlab="Ranking", ylab="PPV", ylim=c(0, 1))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
plot(result.rho1[, 1], result.rho1[, 2], type="l", col="red", xlab="Ranking", ylab="PPV", ylim=c(0, 1), xlim=c(0, 500))
par(new=TRUE)
plot(result.rho1[, 5], result.rho1[, 6], type="l", col="green", xlab="Ranking", ylab="PPV", ylim=c(0, 1), xlim=c(0, 500))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original"), text.col = c("red", "green"), horiz=TRUE)
# Retrieve MT-eQTLs result
MTeQTLs <-read.table(file="MT-eQTLs.txt", header=T)
head(MTeQTLs)
## SNP gene isEQTL.Liver isEQTL.lung marginalP.Liver
## 1 rs6269442 1424963_at 0 0 0.9260306
## 2 rs6365999 1424963_at 0 0 0.9260306
## 3 rs6376963 1424963_at 0 0 0.9306985
## 4 rs3677817 1424963_at 0 0 0.9329117
## 5 rs6269442 1424964_at 0 0 0.9438114
## 6 rs6365999 1424964_at 0 0 0.9438114
## marginalP.lung
## 1 0.9027543
## 2 0.9027543
## 3 0.9104129
## 4 0.9140456
## 5 0.9219186
## 6 0.9219186
mouse430aensembl_id<-read.table(file="2015-12-07 mouse430aensembl_id.txt", header=T)
MTeQTLs<-merge(MTeQTLs, mouse430aensembl_id, by.x = "gene", by.y="probe_id")
MTeQTLs.min <- data.table(MTeQTLs, key=c('ensembl_id', "marginalP.Liver"))
MTeQTLs.min <-MTeQTLs.min[J(unique(ensembl_id)),mult="first"]
merged.eQTL <- merge(new.result.df2.rho1, MTeQTLs.min, by ="ensembl_id")
merged.eQTL$rank.marginalP.Liver <- rank(merged.eQTL$marginalP.Liver)
merged.result <- matrix(, nrow(merged.eQTL), 12)
colnames(merged.result)<-c("bayrank","bayppv","bay_TPR","bay_FPR", "orirank","orippv","ori_TPR","ori_FPR", "MTrank","MTppv","MT_TPR","MT_FPR" )
for (i in 1:nrow(merged.eQTL))
{
newdata1 <- subset(merged.eQTL, rank.p.below.0 <= i)
overlap.newdata1 <- newdata1[newdata1$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
merged.result[i, 1] <- i
merged.result[i, 2] <- nrow(overlap.newdata1)/nrow(newdata1)
merged.result[i, 3] <- nrow(overlap.newdata1)/nrow(liver.ASE.ensembl)
merged.result[i, 4] <- (nrow(newdata1)-nrow(overlap.newdata1)) / (nrow(merged.eQTL)-nrow(liver.ASE.ensembl))
newdata2 <- subset(merged.eQTL, rank.liver.pvalue <= i)
overlap.newdata2 <- newdata2[newdata2$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
merged.result[i, 5] <- i
merged.result[i, 6] <- nrow(overlap.newdata2)/nrow(newdata2)
merged.result[i, 7] <- nrow(overlap.newdata2)/nrow(liver.ASE.ensembl)
merged.result[i, 8] <- (nrow(newdata2)-nrow(overlap.newdata2)) / (nrow(merged.eQTL)-nrow(liver.ASE.ensembl))
newdata3 <- subset(merged.eQTL, rank.marginalP.Liver <= i)
overlap.newdata3 <- newdata3[newdata3$ensembl_id %in% liver.ASE.ensembl$ensembl_gene_id, ]
merged.result[i, 9] <- i
merged.result[i, 10] <- nrow(overlap.newdata3)/nrow(newdata3)
merged.result[i, 11] <- nrow(overlap.newdata3)/nrow(liver.ASE.ensembl)
merged.result[i, 12] <- (nrow(newdata3)-nrow(overlap.newdata3)) /(nrow(merged.eQTL)-nrow(liver.ASE.ensembl))
}
head(merged.result)
## bayrank bayppv bay_TPR bay_FPR orirank orippv
## [1,] 1 NaN 0.000000000 0.0000000000 1 0.0000000
## [2,] 2 0.0000000 0.000000000 0.0002039360 2 0.0000000
## [3,] 3 0.0000000 0.000000000 0.0002039360 3 0.0000000
## [4,] 4 0.3333333 0.005128205 0.0002039360 4 0.3333333
## [5,] 5 0.2500000 0.005128205 0.0003059039 5 0.2500000
## [6,] 6 0.2000000 0.005128205 0.0004078719 6 0.2000000
## ori_TPR ori_FPR MTrank MTppv MT_TPR MT_FPR
## [1,] 0.000000000 0.0001019680 1 0.0000000 0.000000000 0.0001019680
## [2,] 0.000000000 0.0002039360 2 0.5000000 0.005128205 0.0001019680
## [3,] 0.000000000 0.0002039360 3 0.6666667 0.010256410 0.0001019680
## [4,] 0.005128205 0.0002039360 4 0.5000000 0.010256410 0.0002039360
## [5,] 0.005128205 0.0003059039 5 0.6000000 0.015384615 0.0002039360
## [6,] 0.005128205 0.0004078719 6 0.5000000 0.015384615 0.0003059039
tail(merged.result)
## bayrank bayppv bay_TPR bay_FPR orirank orippv
## [9997,] 9997 0.01895785 0.9179487 0.9445294 9997 0.01895785
## [9998,] 9998 0.01895584 0.9179487 0.9446314 9998 0.01895584
## [9999,] 9999 0.01895383 0.9179487 0.9447334 9999 0.01895383
## [10000,] 10000 0.01895183 0.9179487 0.9448353 10000 0.01895183
## [10001,] 10001 0.01894982 0.9179487 0.9449373 10001 0.01894982
## [10002,] 10002 0.01894781 0.9179487 0.9450393 10002 0.01894781
## ori_TPR ori_FPR MTrank MTppv MT_TPR MT_FPR
## [9997,] 0.9179487 0.9445294 9997 0.01840552 0.9435897 1.000612
## [9998,] 0.9179487 0.9446314 9998 0.01840368 0.9435897 1.000714
## [9999,] 0.9179487 0.9447334 9999 0.01840184 0.9435897 1.000816
## [10000,] 0.9179487 0.9448353 10000 0.01840000 0.9435897 1.000918
## [10001,] 0.9179487 0.9449373 10001 0.01839816 0.9435897 1.001020
## [10002,] 0.9179487 0.9450393 10002 0.01839632 0.9435897 1.001122
# ploting "True positive rate"
plot(merged.result[, 1], merged.result[, 3], type="l", col="red", xlab="Ranking", ylab="TPR", ylim=c(0, 1) )
par(new=TRUE)
plot( merged.result[, 1], merged.result[, 7], type="l", col="green", xlab="Ranking", ylab="TPR", ylim=c(0, 1))
par(new=TRUE)
plot(merged.result[, 1], merged.result[, 11], type="l", col="blue", xlab="Ranking", ylab="TPR", ylim=c(0, 1))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original", "MT"), text.col = c("red", "green", "blue"), horiz=TRUE)
plot(merged.result[, 1], merged.result[, 3], type="l", col="red", xlab="Ranking", ylab="TPR", ylim=c(0, 1), xlim=c(0, 300) )
par(new=TRUE)
plot( merged.result[, 1], merged.result[, 7], type="l", col="green", xlab="Ranking", ylab="TPR", ylim=c(0, 1), xlim=c(0, 300))
par(new=TRUE)
plot(merged.result[, 1], merged.result[, 11], type="l", col="blue", xlab="Ranking", ylab="TPR", ylim=c(0, 1), xlim=c(0, 300))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original", "MT"), text.col = c("red", "green", "blue"), horiz=TRUE)
plot(merged.result[, 4], merged.result[, 3], type="l", col="red", xlab="False positive rate", ylab="True positive rate", ylim=c(0, 1) )
par(new=TRUE)
plot( merged.result[, 8], merged.result[, 7], type="l", col="green", xlab="False positive rate", ylab="True positive rate", ylim=c(0, 1))
par(new=TRUE)
plot(merged.result[, 12], merged.result[, 11], type="l", col="blue", xlab="False positive rate", ylab="True positive rate", ylim=c(0, 1))
legend("bottomright", cex = .75, inset=.05, c("Bayesian","Original", "MT"), text.col = c("red", "green", "blue"), horiz=TRUE)
title(main = "ROC curve")
