# Prepare liver gene location and snp location for eQTL analysis
rm(list=ls())
setwd("/Volumes/Transcend/Thesis_project/1.1.liver.gene.snp")
library(moe430a.db)
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
##
## IQR, mad, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, as.vector, cbind,
## colnames, do.call, duplicated, eval, evalq, Filter, Find, get,
## grep, grepl, intersect, is.unsorted, lapply, lengths, Map,
## mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
## setdiff, sort, table, tapply, union, unique, unlist, unsplit
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Loading required package: S4Vectors
## Loading required package: org.Mm.eg.db
## Loading required package: DBI
##
##
x <- moe430aENSEMBL
# Get the entrez gene IDs that are mapped to an Ensembl ID
mapped_genes <- mappedkeys(x)
#For the reverse map ENSEMBL2PROBE:
# Convert to a list
xx <- toTable(moe430aENSEMBL2PROBE)
head(xx)
## probe_id ensembl_id
## 1 1417246_at ENSMUSG00000030359
## 2 1421666_a_at ENSMUSG00000020804
## 3 1416936_at ENSMUSG00000025375
## 4 1421839_at ENSMUSG00000015243
## 5 1421840_at ENSMUSG00000015243
## 6 1450392_at ENSMUSG00000015243
mouse430aensembl_id<-data.frame(xx)
head(mouse430aensembl_id)
## probe_id ensembl_id
## 1 1417246_at ENSMUSG00000030359
## 2 1421666_a_at ENSMUSG00000020804
## 3 1416936_at ENSMUSG00000025375
## 4 1421839_at ENSMUSG00000015243
## 5 1421840_at ENSMUSG00000015243
## 6 1450392_at ENSMUSG00000015243
write.table(mouse430aensembl_id,file="2015-12-07 mouse430aensembl_id.txt", sep="\t", row.names=FALSE, quote=FALSE)
# Get the probe identifiers that are mapped to chromosome locations
xloc <- moe430aCHRLOC
mapped_probesloc <- mappedkeys(xloc)
# Convert to a list
xxloc <- toTable(xloc[mapped_probesloc])
# Per discussion with mentors, take absulute value on end_location
xxloc$start_location<-abs(xxloc$start_location)
head(xxloc)
## probe_id start_location Chromosome
## 1 1415670_at 87887971 6
## 2 1415670_at 87887939 6
## 3 1415671_at 105524469 8
## 4 1415672_at 23241325 8
## 5 1415673_at 129765557 5
## 6 1415674_a_at 44403758 9
xlocend <- moe430aCHRLOCEND
# Get the probe identifiers that are mapped to chromosome locations
mapped_probes <- mappedkeys(xlocend)
# Convert to a table
xxlocend <- toTable(xlocend[mapped_probes])
head(xxlocend)
## probe_id end_location Chromosome
## 1 1415670_at 87890759 6
## 2 1415670_at 87913594 6
## 3 1415671_at -105566040 8
## 4 1415672_at -23257080 8
## 5 1415673_at -129787253 5
## 6 1415674_a_at -44407548 9
# Per discussion with mentors, take absulute value on end_location
xxlocend$end_location<-abs(xxlocend$end_location)
dim(xxlocend)
## [1] 26903 3
# combine gene location information to one table
mouse430achrloc<-cbind(xxloc, xxlocend$end_location)
mouse430achrloc <- mouse430achrloc[c("probe_id", "Chromosome", "start_location", "xxlocend$end_location")]
head(mouse430achrloc)
## probe_id Chromosome start_location xxlocend$end_location
## 1 1415670_at 6 87887971 87890759
## 2 1415670_at 6 87887939 87913594
## 3 1415671_at 8 105524469 105566040
## 4 1415672_at 8 23241325 23257080
## 5 1415673_at 5 129765557 129787253
## 6 1415674_a_at 9 44403758 44407548
library(plyr)
##
## Attaching package: 'plyr'
## The following object is masked from 'package:IRanges':
##
## desc
## The following object is masked from 'package:S4Vectors':
##
## rename
mouse430achrloc<-rename(mouse430achrloc, c("xxlocend$end_location"="end_location"))
# delete duplicated probe_id (difference of )
mouse430achrloc<-mouse430achrloc[!duplicated(mouse430achrloc$probe_id),]
# sort by probe_id
mouse430achrloc <-mouse430achrloc[order(mouse430achrloc$probe_id),]
head(mouse430achrloc)
## probe_id Chromosome start_location end_location
## 1 1415670_at 6 87887971 87890759
## 3 1415671_at 8 105524469 105566040
## 4 1415672_at 8 23241325 23257080
## 5 1415673_at 5 129765557 129787253
## 6 1415674_a_at 9 44403758 44407548
## 7 1415675_at 2 32570857 32573571
# write gene location for matrix eqtl analysis
liver.gene.loc <- mouse430achrloc
write.table(liver.gene.loc,file="2015-12-07 liver.gene.loc.txt", sep="\t", row.names=FALSE, quote=FALSE)
# Prepare SNP location for liver eQTL anlaysis
BXD.geno<-read.table(file="BXD-3.geno.txt", header=T)
head(BXD.geno)
## Chr Locus cM Mb BXD1 BXD2 BXD5 BXD6 BXD8 BXD9 BXD11 BXD12
## 1 1 rs6269442 0.000 3.482275 B B D D D B B D
## 2 1 rs6365999 0.000 4.811062 B B D D D B B D
## 3 1 rs6376963 0.895 5.008089 B B D D D B B D
## 4 1 rs3677817 1.185 5.176058 B B D D D B B D
## 5 1 rs8236463 2.081 5.579193 B B D D D B D D
## 6 1 rs6333200 2.081 6.217921 B B D D D B D D
## BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD20 BXD21 BXD22 BXD23 BXD24 BXD24a
## 1 B B D D B D D D D B B B
## 2 B B D D B D D D D B B B
## 3 B B D D B D D D D B B B
## 4 B B D D B D D D D B B B
## 5 B B D D B D D D D B B B
## 6 B B D D B D D D D B B B
## BXD25 BXD27 BXD28 BXD29 BXD30 BXD31 BXD32 BXD33 BXD34 BXD35 BXD36 BXD37
## 1 D B D D B B B B B B B B
## 2 D B D D B B B B B B B B
## 3 D B D D B B B B B B B B
## 4 D B D D B B B B B B B B
## 5 D B D D B B B B B B B B
## 6 D B D D B B B B B B B B
## BXD38 BXD39 BXD40 BXD41 BXD42 BXD43 BXD44 BXD45 BXD48 BXD48a BXD49 BXD50
## 1 B D B D B B D B B B H H
## 2 B D B D B B D B B B H H
## 3 B D B D B D D B B B H H
## 4 B D B D B D D B B B H H
## 5 B D B D B D D B B B H H
## 6 B D B D B D D B B B H H
## BXD51 BXD52 BXD53 BXD54 BXD55 BXD56 BXD59 BXD60 BXD61 BXD62 BXD63 BXD64
## 1 B D B B H H B B D D D D
## 2 B D B B H H B B D D D D
## 3 B B B B H H B B D D D D
## 4 B B B B H H B B D D D D
## 5 B B B B H H B B D D D D
## 6 B B B B H H B B D D D D
## BXD65 BXD65a BXD65b BXD66 BXD67 BXD68 BXD69 BXD70 BXD71 BXD72 BXD73
## 1 D B D B B H B B B B D
## 2 D B D B B H B B B B D
## 3 B B D B B H B B B B D
## 4 B B D B B D B B B B D
## 5 B B D B B D B B B B D
## 6 B B D B B D B B B B D
## BXD73a BXD73b BXD74 BXD75 BXD76 BXD77 BXD78 BXD79 BXD81 BXD83 BXD84
## 1 D D B D B D D D D H B
## 2 D D B D B D D D D H B
## 3 D D B D B D D D D H B
## 4 D D B D B D D B D H B
## 5 D D B D B B D B B H B
## 6 D D B D B B D B B H B
## BXD85 BXD86 BXD87 BXD88 BXD89 BXD90 BXD91 BXD93 BXD94 BXD95 BXD98 BXD99
## 1 D D B D B B D B D D B B
## 2 D D B D B B D B D D B B
## 3 D D B D B B D B D D B B
## 4 D D B D B B D B D D B B
## 5 D D B D B B D B D D B B
## 6 D D B D B B D B D D B B
## BXD100 BXD101 BXD102
## 1 B B B
## 2 B B U
## 3 B B U
## 4 B B U
## 5 B B U
## 6 B B U
dim(BXD.geno)
## [1] 3811 97
# recode SNP, "B to 0, D to 1"
BXD.geno1 <- as.data.frame(sapply(BXD.geno,gsub,pattern="B",replacement="0"))
# BXD.geno2 <- as.data.frame(sapply(BXD.geno1,gsub,pattern="H",replacement="1"))
# BXD.geno3 <- as.data.frame(sapply(BXD.geno2,gsub,pattern="D",replacement="2"))
BXD.geno2 <- as.data.frame(sapply(BXD.geno1,gsub,pattern="H",replacement="NA"))
BXD.geno3 <- as.data.frame(sapply(BXD.geno2,gsub,pattern="D",replacement="1"))
BXD.geno4 <- as.data.frame(sapply(BXD.geno3,gsub,pattern="U",replacement="NA"))
head(BXD.geno4)
## Chr Locus cM Mb BXD1 BXD2 BXD5 BXD6 BXD8 BXD9 BXD11 BXD12
## 1 1 rs6269442 0 3.482275 0 0 1 1 1 0 0 1
## 2 1 rs6365999 0 4.811062 0 0 1 1 1 0 0 1
## 3 1 rs6376963 0.895 5.008089 0 0 1 1 1 0 0 1
## 4 1 rs3677817 1.185 5.176058 0 0 1 1 1 0 0 1
## 5 1 rs8236463 2.081 5.579193 0 0 1 1 1 0 1 1
## 6 1 rs6333200 2.081 6.217921 0 0 1 1 1 0 1 1
## BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD20 BXD21 BXD22 BXD23 BXD24 BXD24a
## 1 0 0 1 1 0 1 1 1 1 0 0 0
## 2 0 0 1 1 0 1 1 1 1 0 0 0
## 3 0 0 1 1 0 1 1 1 1 0 0 0
## 4 0 0 1 1 0 1 1 1 1 0 0 0
## 5 0 0 1 1 0 1 1 1 1 0 0 0
## 6 0 0 1 1 0 1 1 1 1 0 0 0
## BXD25 BXD27 BXD28 BXD29 BXD30 BXD31 BXD32 BXD33 BXD34 BXD35 BXD36 BXD37
## 1 1 0 1 1 0 0 0 0 0 0 0 0
## 2 1 0 1 1 0 0 0 0 0 0 0 0
## 3 1 0 1 1 0 0 0 0 0 0 0 0
## 4 1 0 1 1 0 0 0 0 0 0 0 0
## 5 1 0 1 1 0 0 0 0 0 0 0 0
## 6 1 0 1 1 0 0 0 0 0 0 0 0
## BXD38 BXD39 BXD40 BXD41 BXD42 BXD43 BXD44 BXD45 BXD48 BXD48a BXD49 BXD50
## 1 0 1 0 1 0 0 1 0 0 0 NA NA
## 2 0 1 0 1 0 0 1 0 0 0 NA NA
## 3 0 1 0 1 0 1 1 0 0 0 NA NA
## 4 0 1 0 1 0 1 1 0 0 0 NA NA
## 5 0 1 0 1 0 1 1 0 0 0 NA NA
## 6 0 1 0 1 0 1 1 0 0 0 NA NA
## BXD51 BXD52 BXD53 BXD54 BXD55 BXD56 BXD59 BXD60 BXD61 BXD62 BXD63 BXD64
## 1 0 1 0 0 NA NA 0 0 1 1 1 1
## 2 0 1 0 0 NA NA 0 0 1 1 1 1
## 3 0 0 0 0 NA NA 0 0 1 1 1 1
## 4 0 0 0 0 NA NA 0 0 1 1 1 1
## 5 0 0 0 0 NA NA 0 0 1 1 1 1
## 6 0 0 0 0 NA NA 0 0 1 1 1 1
## BXD65 BXD65a BXD65b BXD66 BXD67 BXD68 BXD69 BXD70 BXD71 BXD72 BXD73
## 1 1 0 1 0 0 NA 0 0 0 0 1
## 2 1 0 1 0 0 NA 0 0 0 0 1
## 3 0 0 1 0 0 NA 0 0 0 0 1
## 4 0 0 1 0 0 1 0 0 0 0 1
## 5 0 0 1 0 0 1 0 0 0 0 1
## 6 0 0 1 0 0 1 0 0 0 0 1
## BXD73a BXD73b BXD74 BXD75 BXD76 BXD77 BXD78 BXD79 BXD81 BXD83 BXD84
## 1 1 1 0 1 0 1 1 1 1 NA 0
## 2 1 1 0 1 0 1 1 1 1 NA 0
## 3 1 1 0 1 0 1 1 1 1 NA 0
## 4 1 1 0 1 0 1 1 0 1 NA 0
## 5 1 1 0 1 0 0 1 0 0 NA 0
## 6 1 1 0 1 0 0 1 0 0 NA 0
## BXD85 BXD86 BXD87 BXD88 BXD89 BXD90 BXD91 BXD93 BXD94 BXD95 BXD98 BXD99
## 1 1 1 0 1 0 0 1 0 1 1 0 0
## 2 1 1 0 1 0 0 1 0 1 1 0 0
## 3 1 1 0 1 0 0 1 0 1 1 0 0
## 4 1 1 0 1 0 0 1 0 1 1 0 0
## 5 1 1 0 1 0 0 1 0 1 1 0 0
## 6 1 1 0 1 0 0 1 0 1 1 0 0
## BXD100 BXD101 BXD102
## 1 0 0 0
## 2 0 0 NA
## 3 0 0 NA
## 4 0 0 NA
## 5 0 0 NA
## 6 0 0 NA
BXD.geno.loc<- BXD.geno4[, c(2, 1, 4)]
# Transform SNP position from Mb to bp
BXD.geno.loc$Mb<-as.numeric(as.character(BXD.geno.loc$Mb))
BXD.geno.loc[, 3]<-BXD.geno.loc[, 3] * 1000000
library(plyr)
BXD.geno.loc<-rename(BXD.geno.loc, c("Mb"="position"))
head(BXD.geno.loc)
## Locus Chr position
## 1 rs6269442 1 3482275
## 2 rs6365999 1 4811062
## 3 rs6376963 1 5008089
## 4 rs3677817 1 5176058
## 5 rs8236463 1 5579193
## 6 rs6333200 1 6217921
BXD.geno.SNP<-BXD.geno4[, c(2, 5:97)]
head(BXD.geno.SNP)
## Locus BXD1 BXD2 BXD5 BXD6 BXD8 BXD9 BXD11 BXD12 BXD13 BXD14 BXD15
## 1 rs6269442 0 0 1 1 1 0 0 1 0 0 1
## 2 rs6365999 0 0 1 1 1 0 0 1 0 0 1
## 3 rs6376963 0 0 1 1 1 0 0 1 0 0 1
## 4 rs3677817 0 0 1 1 1 0 0 1 0 0 1
## 5 rs8236463 0 0 1 1 1 0 1 1 0 0 1
## 6 rs6333200 0 0 1 1 1 0 1 1 0 0 1
## BXD16 BXD18 BXD19 BXD20 BXD21 BXD22 BXD23 BXD24 BXD24a BXD25 BXD27 BXD28
## 1 1 0 1 1 1 1 0 0 0 1 0 1
## 2 1 0 1 1 1 1 0 0 0 1 0 1
## 3 1 0 1 1 1 1 0 0 0 1 0 1
## 4 1 0 1 1 1 1 0 0 0 1 0 1
## 5 1 0 1 1 1 1 0 0 0 1 0 1
## 6 1 0 1 1 1 1 0 0 0 1 0 1
## BXD29 BXD30 BXD31 BXD32 BXD33 BXD34 BXD35 BXD36 BXD37 BXD38 BXD39 BXD40
## 1 1 0 0 0 0 0 0 0 0 0 1 0
## 2 1 0 0 0 0 0 0 0 0 0 1 0
## 3 1 0 0 0 0 0 0 0 0 0 1 0
## 4 1 0 0 0 0 0 0 0 0 0 1 0
## 5 1 0 0 0 0 0 0 0 0 0 1 0
## 6 1 0 0 0 0 0 0 0 0 0 1 0
## BXD41 BXD42 BXD43 BXD44 BXD45 BXD48 BXD48a BXD49 BXD50 BXD51 BXD52 BXD53
## 1 1 0 0 1 0 0 0 NA NA 0 1 0
## 2 1 0 0 1 0 0 0 NA NA 0 1 0
## 3 1 0 1 1 0 0 0 NA NA 0 0 0
## 4 1 0 1 1 0 0 0 NA NA 0 0 0
## 5 1 0 1 1 0 0 0 NA NA 0 0 0
## 6 1 0 1 1 0 0 0 NA NA 0 0 0
## BXD54 BXD55 BXD56 BXD59 BXD60 BXD61 BXD62 BXD63 BXD64 BXD65 BXD65a
## 1 0 NA NA 0 0 1 1 1 1 1 0
## 2 0 NA NA 0 0 1 1 1 1 1 0
## 3 0 NA NA 0 0 1 1 1 1 0 0
## 4 0 NA NA 0 0 1 1 1 1 0 0
## 5 0 NA NA 0 0 1 1 1 1 0 0
## 6 0 NA NA 0 0 1 1 1 1 0 0
## BXD65b BXD66 BXD67 BXD68 BXD69 BXD70 BXD71 BXD72 BXD73 BXD73a BXD73b
## 1 1 0 0 NA 0 0 0 0 1 1 1
## 2 1 0 0 NA 0 0 0 0 1 1 1
## 3 1 0 0 NA 0 0 0 0 1 1 1
## 4 1 0 0 1 0 0 0 0 1 1 1
## 5 1 0 0 1 0 0 0 0 1 1 1
## 6 1 0 0 1 0 0 0 0 1 1 1
## BXD74 BXD75 BXD76 BXD77 BXD78 BXD79 BXD81 BXD83 BXD84 BXD85 BXD86 BXD87
## 1 0 1 0 1 1 1 1 NA 0 1 1 0
## 2 0 1 0 1 1 1 1 NA 0 1 1 0
## 3 0 1 0 1 1 1 1 NA 0 1 1 0
## 4 0 1 0 1 1 0 1 NA 0 1 1 0
## 5 0 1 0 0 1 0 0 NA 0 1 1 0
## 6 0 1 0 0 1 0 0 NA 0 1 1 0
## BXD88 BXD89 BXD90 BXD91 BXD93 BXD94 BXD95 BXD98 BXD99 BXD100 BXD101
## 1 1 0 0 1 0 1 1 0 0 0 0
## 2 1 0 0 1 0 1 1 0 0 0 0
## 3 1 0 0 1 0 1 1 0 0 0 0
## 4 1 0 0 1 0 1 1 0 0 0 0
## 5 1 0 0 1 0 1 1 0 0 0 0
## 6 1 0 0 1 0 1 1 0 0 0 0
## BXD102
## 1 0
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
# Reformat mouse liver gene expression for Matrix eqtl analysis
mouse.liver.expression<-read.table("GN373_GSE16780_UCLA_Hybrid_MDP_Liver_Affy_HT_M430A_Sep11_RMA_Z-Score_Average.txt", comment.char ="#", header = TRUE, sep = "\t", )
head(mouse.liver.expression)
## ProbeSet Gene.Symbol C57BL.6J DBA.2J BXD1 BXD2 BXD5 BXD6
## 1 1422417_at_A X99300 6.320 6.283 6.280 6.276 6.487 6.392
## 2 1426607_at_A 3110070M22Rik 7.240 6.756 6.567 7.219 6.718 6.542
## 3 1434698_at_A Adamts17 10.380 10.422 10.653 10.430 10.871 10.544
## 4 1435792_at_A Csprs 6.809 7.670 6.980 6.806 6.898 6.838
## 5 1437243_at_A D2hgdh 7.876 8.139 8.113 7.938 8.460 8.249
## 6 1438225_x_at_A Tram1 6.867 6.842 6.846 6.970 6.989 7.007
## BXD8 BXD9 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18 BXD19
## 1 6.276 6.358 6.331 6.398 6.294 6.363 6.309 6.261 6.174 6.303
## 2 6.549 6.875 6.584 6.585 7.098 7.087 7.085 6.586 7.903 7.232
## 3 10.300 10.522 10.580 10.745 10.599 10.741 10.764 10.697 10.653 10.612
## 4 7.081 6.894 6.731 7.302 6.955 7.315 6.812 7.054 6.958 7.004
## 5 7.742 8.066 7.872 8.115 8.348 8.036 7.954 8.087 7.978 8.293
## 6 6.838 7.054 6.811 6.924 7.035 6.977 6.800 6.879 7.109 6.914
## BXD20 BXD21 BXD24 BXD27 BXD28 BXD29 BXD31 BXD32 BXD33 BXD34
## 1 6.339 6.256 6.306 6.236 6.351 6.425 6.346 6.351 6.402 6.342
## 2 6.594 7.417 6.645 7.209 7.756 7.169 6.465 6.612 6.642 7.407
## 3 10.493 10.785 10.664 10.761 10.536 10.619 10.855 10.765 10.475 10.738
## 4 7.145 7.543 6.860 6.921 6.874 7.143 7.043 7.075 7.147 6.776
## 5 8.049 7.898 8.145 7.598 8.072 7.942 8.104 8.368 7.697 8.072
## 6 6.821 6.879 6.957 6.704 7.082 6.786 6.887 7.039 6.937 6.999
## BXD36 BXD38 BXD39 BXD40 BXD42 A.J AXB1 AXB2 AXB4 AXB5
## 1 6.240 6.289 6.244 6.396 6.337 6.382 6.337 6.338 6.475 6.294
## 2 6.582 7.261 7.400 7.932 6.412 7.920 7.339 8.879 7.618 8.515
## 3 10.754 10.774 10.643 10.619 10.385 10.642 10.741 10.671 10.702 10.712
## 4 7.340 7.459 7.111 7.126 7.053 6.801 6.993 6.970 6.804 6.666
## 5 7.864 8.020 8.470 7.978 8.226 8.170 8.008 8.055 8.268 7.806
## 6 6.922 6.689 6.832 7.113 6.801 6.820 6.930 6.798 6.886 6.725
## AXB6 AXB8 AXB10 AXB12 AXB15 AXB23 AXB24 BXA1 BXA2 BXA4
## 1 6.187 6.312 6.408 6.377 6.233 6.407 6.264 6.191 6.357 6.251
## 2 7.581 7.543 7.227 8.626 7.315 8.915 7.142 7.352 7.735 7.546
## 3 10.654 10.712 10.806 10.690 10.761 10.928 10.735 10.360 10.453 10.820
## 4 6.985 6.768 6.950 6.865 6.755 6.661 6.576 6.668 6.806 6.800
## 5 7.918 8.134 8.556 8.360 7.963 8.147 8.136 8.045 8.477 8.264
## 6 6.808 6.826 7.049 6.826 6.863 6.880 6.948 7.020 6.853 6.804
## BXA7 BXA11 BXA12 BXA13 BXA14 BXA16 BXA24 BXA25 BXA26 AKR.J
## 1 6.400 6.509 6.298 6.313 6.259 6.290 6.197 6.394 6.308 6.308
## 2 8.422 9.367 7.175 9.066 7.387 7.322 8.070 7.818 9.088 8.320
## 3 10.741 10.654 10.560 10.837 10.677 10.658 10.609 10.552 10.662 10.532
## 4 6.786 6.995 6.965 6.732 6.825 6.848 6.643 6.742 6.914 7.927
## 5 8.191 8.258 8.192 8.106 8.234 8.265 8.072 8.059 8.067 8.302
## 6 6.911 6.845 6.862 6.877 7.082 6.901 6.952 6.880 6.831 6.996
## C3H.HeJ BXH2 BXH4 BXH6 BXH7 BXH8 BXH9 BXH10 BXH14 BXH19
## 1 6.361 6.284 6.363 6.390 6.257 6.302 6.396 6.366 6.382 6.310
## 2 6.662 6.442 7.832 7.432 6.621 6.393 7.497 6.325 6.657 6.428
## 3 10.798 10.229 10.537 10.588 10.707 10.558 10.473 10.510 10.634 10.775
## 4 7.060 6.622 7.069 6.788 6.892 7.065 6.875 6.803 7.108 7.087
## 5 8.267 7.501 8.116 8.241 8.216 8.289 8.553 8.268 8.393 8.415
## 6 6.945 6.814 6.832 6.944 6.932 6.892 6.966 6.897 6.944 7.133
## BXH20 BXH21 B6cC3.1 BALB.cJ LG.J NOD.ShiLtJ AXB13.14 AXB18.19.20
## 1 6.532 6.329 6.300 6.313 6.328 6.312 6.246 6.314
## 2 6.457 6.551 6.616 8.518 7.353 8.062 7.515 7.359
## 3 11.015 10.762 10.493 10.651 10.560 10.792 10.725 10.829
## 4 6.841 6.900 7.019 7.028 7.591 7.648 6.977 6.841
## 5 7.568 8.474 8.177 8.507 8.187 8.348 8.189 8.276
## 6 6.629 6.926 6.961 6.839 6.823 6.793 6.814 6.857
## BXA8.17 X129X1.SvJ BTBRT...tf.J BUB.BnJ C57L.J C58.J CAST.EiJ CBA.J
## 1 6.343 6.271 6.381 6.371 6.302 6.202 6.335 6.346
## 2 7.419 6.874 6.552 7.808 7.484 7.269 6.659 6.443
## 3 10.762 10.520 10.716 10.780 10.580 10.434 10.005 10.692
## 4 6.900 6.744 6.816 7.116 6.840 6.650 6.648 7.640
## 5 8.189 8.085 7.873 8.060 8.126 8.084 7.888 8.425
## 6 6.858 7.036 6.924 6.964 6.998 6.833 6.754 6.948
## CE.J FVB.NJ I.LnJ KK.HlJ LP.J MA.MyJ NON.LtJ NZB.BlNJ NZW.LacJ
## 1 6.255 6.312 6.352 6.450 6.242 6.274 6.209 6.296 6.302
## 2 6.371 8.073 7.578 7.830 6.501 6.752 8.257 7.552 7.472
## 3 10.723 10.734 10.308 10.539 10.493 10.572 10.463 10.531 10.514
## 4 6.969 6.814 6.741 7.268 6.631 6.766 6.988 6.848 6.835
## 5 8.111 8.681 8.054 8.699 7.854 7.807 8.133 8.028 8.108
## 6 6.899 7.008 7.022 6.854 6.861 6.842 6.961 7.065 7.026
## PL.J RIIIS.J SEA.GnJ SM.J SWR.J BXD24a X
## 1 6.326 6.438 6.288 6.248 6.321 6.317 NA
## 2 9.147 6.566 7.410 6.545 8.019 6.817 NA
## 3 10.487 10.591 10.478 10.607 10.465 10.654 NA
## 4 7.252 6.809 7.094 7.320 6.905 7.001 NA
## 5 7.149 7.941 8.259 7.934 7.686 8.211 NA
## 6 6.796 6.908 6.958 6.806 6.822 6.793 NA
mouse.liver.expression<- as.data.frame(sapply(mouse.liver.expression,gsub,pattern="_at_A",replacement="_at"))
# Creat the strain library with known SNP
BXD.geno.SNP.library<-colnames(BXD.geno.SNP)
length(BXD.geno.SNP.library)
## [1] 94
mouse.liver.expression.eqtl<-mouse.liver.expression[, which(colnames(mouse.liver.expression) %in% BXD.geno.SNP.library)]
# reorder stain column names
mouse.liver.expression.eqtl<-mouse.liver.expression.eqtl[,order(colnames(mouse.liver.expression.eqtl))]
mouse.liver.expression.eqtl$ProbeSet<-mouse.liver.expression$ProbeSet
mouse.liver.expression.eqtl<-mouse.liver.expression.eqtl[, c(31, 1:30)]
head(mouse.liver.expression.eqtl)
## ProbeSet BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18
## 1 1422417_at 6.28 6.331 6.398 6.294 6.363 6.309 6.261 6.174
## 2 1426607_at 6.567 6.584 6.585 7.098 7.087 7.085 6.586 7.903
## 3 1434698_at 10.653 10.58 10.745 10.599 10.741 10.764 10.697 10.653
## 4 1435792_at 6.98 6.731 7.302 6.955 7.315 6.812 7.054 6.958
## 5 1437243_at 8.113 7.872 8.115 8.348 8.036 7.954 8.087 7.978
## 6 1438225_x_at 6.846 6.811 6.924 7.035 6.977 6.8 6.879 7.109
## BXD19 BXD2 BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31
## 1 6.303 6.276 6.339 6.256 6.306 6.317 6.236 6.351 6.425 6.346
## 2 7.232 7.219 6.594 7.417 6.645 6.817 7.209 7.756 7.169 6.465
## 3 10.612 10.43 10.493 10.785 10.664 10.654 10.761 10.536 10.619 10.855
## 4 7.004 6.806 7.145 7.543 6.86 7.001 6.921 6.874 7.143 7.043
## 5 8.293 7.938 8.049 7.898 8.145 8.211 7.598 8.072 7.942 8.104
## 6 6.914 6.97 6.821 6.879 6.957 6.793 6.704 7.082 6.786 6.887
## BXD32 BXD33 BXD34 BXD36 BXD38 BXD39 BXD40 BXD42 BXD5 BXD6
## 1 6.351 6.402 6.342 6.24 6.289 6.244 6.396 6.337 6.487 6.392
## 2 6.612 6.642 7.407 6.582 7.261 7.4 7.932 6.412 6.718 6.542
## 3 10.765 10.475 10.738 10.754 10.774 10.643 10.619 10.385 10.871 10.544
## 4 7.075 7.147 6.776 7.34 7.459 7.111 7.126 7.053 6.898 6.838
## 5 8.368 7.697 8.072 7.864 8.02 8.47 7.978 8.226 8.46 8.249
## 6 7.039 6.937 6.999 6.922 6.689 6.832 7.113 6.801 6.989 7.007
## BXD8 BXD9
## 1 6.276 6.358
## 2 6.549 6.875
## 3 10.3 10.522
## 4 7.081 6.894
## 5 7.742 8.066
## 6 6.838 7.054
# creat strain library with liver expression data
BXD.geno.strain.library<-colnames(mouse.liver.expression.eqtl)
# select SNP on the strains which has gene expression data available
BXD.geno5<-BXD.geno4[, which(colnames(BXD.geno4) %in% BXD.geno.strain.library)]
BXD.geno5<-BXD.geno5[,order(colnames(BXD.geno5))]
BXD.geno5$Locus<-BXD.geno4$Locus
BXD.geno.SNP.eqtl<-BXD.geno5[, c(31, 1:30)]
head(BXD.geno.SNP.eqtl)
## Locus BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD2
## 1 rs6269442 0 0 1 0 0 1 1 0 1 0
## 2 rs6365999 0 0 1 0 0 1 1 0 1 0
## 3 rs6376963 0 0 1 0 0 1 1 0 1 0
## 4 rs3677817 0 0 1 0 0 1 1 0 1 0
## 5 rs8236463 0 1 1 0 0 1 1 0 1 0
## 6 rs6333200 0 1 1 0 0 1 1 0 1 0
## BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31 BXD32 BXD33 BXD34 BXD36
## 1 1 1 0 0 0 1 1 0 0 0 0 0
## 2 1 1 0 0 0 1 1 0 0 0 0 0
## 3 1 1 0 0 0 1 1 0 0 0 0 0
## 4 1 1 0 0 0 1 1 0 0 0 0 0
## 5 1 1 0 0 0 1 1 0 0 0 0 0
## 6 1 1 0 0 0 1 1 0 0 0 0 0
## BXD38 BXD39 BXD40 BXD42 BXD5 BXD6 BXD8 BXD9
## 1 0 1 0 0 1 1 1 0
## 2 0 1 0 0 1 1 1 0
## 3 0 1 0 0 1 1 1 0
## 4 0 1 0 0 1 1 1 0
## 5 0 1 0 0 1 1 1 0
## 6 0 1 0 0 1 1 1 0
dim(BXD.geno.SNP.eqtl)
## [1] 3811 31
# write BXD SNP genotypes for eqtl analysis
write.table(BXD.geno.SNP.eqtl,file="2015-12-07 BXD.geno.SNP.eqtl.for.liver.txt", sep="\t", row.names=FALSE, quote=FALSE)
# check dimensions to make sure they match
dim(BXD.geno.loc)
## [1] 3811 3
head(BXD.geno.loc)
## Locus Chr position
## 1 rs6269442 1 3482275
## 2 rs6365999 1 4811062
## 3 rs6376963 1 5008089
## 4 rs3677817 1 5176058
## 5 rs8236463 1 5579193
## 6 rs6333200 1 6217921
head(BXD.geno.SNP.eqtl)
## Locus BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD2
## 1 rs6269442 0 0 1 0 0 1 1 0 1 0
## 2 rs6365999 0 0 1 0 0 1 1 0 1 0
## 3 rs6376963 0 0 1 0 0 1 1 0 1 0
## 4 rs3677817 0 0 1 0 0 1 1 0 1 0
## 5 rs8236463 0 1 1 0 0 1 1 0 1 0
## 6 rs6333200 0 1 1 0 0 1 1 0 1 0
## BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31 BXD32 BXD33 BXD34 BXD36
## 1 1 1 0 0 0 1 1 0 0 0 0 0
## 2 1 1 0 0 0 1 1 0 0 0 0 0
## 3 1 1 0 0 0 1 1 0 0 0 0 0
## 4 1 1 0 0 0 1 1 0 0 0 0 0
## 5 1 1 0 0 0 1 1 0 0 0 0 0
## 6 1 1 0 0 0 1 1 0 0 0 0 0
## BXD38 BXD39 BXD40 BXD42 BXD5 BXD6 BXD8 BXD9
## 1 0 1 0 0 1 1 1 0
## 2 0 1 0 0 1 1 1 0
## 3 0 1 0 0 1 1 1 0
## 4 0 1 0 0 1 1 1 0
## 5 0 1 0 0 1 1 1 0
## 6 0 1 0 0 1 1 1 0
head(mouse.liver.expression.eqtl)
## ProbeSet BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18
## 1 1422417_at 6.28 6.331 6.398 6.294 6.363 6.309 6.261 6.174
## 2 1426607_at 6.567 6.584 6.585 7.098 7.087 7.085 6.586 7.903
## 3 1434698_at 10.653 10.58 10.745 10.599 10.741 10.764 10.697 10.653
## 4 1435792_at 6.98 6.731 7.302 6.955 7.315 6.812 7.054 6.958
## 5 1437243_at 8.113 7.872 8.115 8.348 8.036 7.954 8.087 7.978
## 6 1438225_x_at 6.846 6.811 6.924 7.035 6.977 6.8 6.879 7.109
## BXD19 BXD2 BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31
## 1 6.303 6.276 6.339 6.256 6.306 6.317 6.236 6.351 6.425 6.346
## 2 7.232 7.219 6.594 7.417 6.645 6.817 7.209 7.756 7.169 6.465
## 3 10.612 10.43 10.493 10.785 10.664 10.654 10.761 10.536 10.619 10.855
## 4 7.004 6.806 7.145 7.543 6.86 7.001 6.921 6.874 7.143 7.043
## 5 8.293 7.938 8.049 7.898 8.145 8.211 7.598 8.072 7.942 8.104
## 6 6.914 6.97 6.821 6.879 6.957 6.793 6.704 7.082 6.786 6.887
## BXD32 BXD33 BXD34 BXD36 BXD38 BXD39 BXD40 BXD42 BXD5 BXD6
## 1 6.351 6.402 6.342 6.24 6.289 6.244 6.396 6.337 6.487 6.392
## 2 6.612 6.642 7.407 6.582 7.261 7.4 7.932 6.412 6.718 6.542
## 3 10.765 10.475 10.738 10.754 10.774 10.643 10.619 10.385 10.871 10.544
## 4 7.075 7.147 6.776 7.34 7.459 7.111 7.126 7.053 6.898 6.838
## 5 8.368 7.697 8.072 7.864 8.02 8.47 7.978 8.226 8.46 8.249
## 6 7.039 6.937 6.999 6.922 6.689 6.832 7.113 6.801 6.989 7.007
## BXD8 BXD9
## 1 6.276 6.358
## 2 6.549 6.875
## 3 10.3 10.522
## 4 7.081 6.894
## 5 7.742 8.066
## 6 6.838 7.054
# write BXD SNP location for eqtl analysis
write.table(BXD.geno.loc,file="2015-12-07 BXD.geno.loc.eqtl.for.liver.txt", sep="\t", row.names=FALSE, quote=FALSE)
liver.probeset.position.library<-mouse430achrloc$probe_id
length(liver.probeset.position.library)
## [1] 20855
# subset mouse liver expression data with known gene location
mouse.liver.expression.eqtl.position<-mouse.liver.expression.eqtl[mouse.liver.expression.eqtl$ProbeSet %in% liver.probeset.position.library,]
mouse.liver.expression.eqtl.position <-mouse.liver.expression.eqtl.position[order(mouse.liver.expression.eqtl.position$ProbeSet),]
head(mouse.liver.expression.eqtl.position)
## ProbeSet BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18
## 18098 1415670_at 10.09 10.2 10.3 10.208 9.83 10.238 9.914 10.348
## 21104 1415671_at 10.932 11.088 11.007 11.02 10.955 11.12 11.012 11.123
## 16745 1415672_at 11.432 11.417 11.442 11.555 11.561 11.318 11.461 11.561
## 17051 1415673_at 7.535 7.382 7.566 7.162 7.403 7.342 7.213 7.581
## 21715 1415674_a_at 9.757 9.972 9.269 9.873 9.354 9.918 9.459 9.655
## 11364 1415675_at 9.029 9.009 9.245 9.282 9.415 9.098 9.06 8.937
## BXD19 BXD2 BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29
## 18098 9.939 9.871 10.077 10.159 9.746 9.89 10.286 10.177 9.959
## 21104 10.922 10.802 10.988 10.969 11.104 10.979 10.905 11.013 11.071
## 16745 11.575 11.426 11.367 11.328 11.499 11.531 11.666 11.5 11.502
## 17051 7.551 7.368 7.251 7.373 7.408 7.473 7.238 7.424 7.475
## 21715 9.544 9.557 9.46 9.322 9.795 9.758 9.628 9.263 9.671
## 11364 8.995 9.104 9.123 9.086 9.087 9.04 8.844 9.173 8.951
## BXD31 BXD32 BXD33 BXD34 BXD36 BXD38 BXD39 BXD40 BXD42
## 18098 9.882 10.059 10.102 10.174 10.022 10.364 9.745 10.074 9.961
## 21104 10.996 10.854 11.084 11.059 10.923 11.053 11.03 11.067 10.891
## 16745 11.446 11.548 11.511 11.49 11.545 11.496 11.516 11.543 11.457
## 17051 7.491 7.419 7.386 7.496 7.369 7.265 7.13 7.149 7.383
## 21715 9.645 9.368 9.741 9.721 9.684 9.387 9.649 9.693 9.793
## 11364 9.131 9.108 9.122 9.153 9.008 9.138 9.126 9.233 9.06
## BXD5 BXD6 BXD8 BXD9
## 18098 10.16 10.069 9.956 10.142
## 21104 10.878 11.043 11.206 10.982
## 16745 11.374 11.504 11.531 11.7
## 17051 7.439 7.233 7.397 7.39
## 21715 9.84 9.519 9.34 9.16
## 11364 9.134 9.081 8.926 9.274
dim(mouse.liver.expression.eqtl.position)
## [1] 20855 31
# write mouse liver gene expression data for eqtl analsis
write.table(mouse.liver.expression.eqtl.position,file="2015-12-07 mouse.liver.expression.eqtl.txt", sep="\t", row.names=FALSE, quote=FALSE)
dim(mouse.liver.expression.eqtl.position)
## [1] 20855 31
dim(liver.gene.loc)
## [1] 20855 4
head(mouse.liver.expression.eqtl.position)
## ProbeSet BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18
## 18098 1415670_at 10.09 10.2 10.3 10.208 9.83 10.238 9.914 10.348
## 21104 1415671_at 10.932 11.088 11.007 11.02 10.955 11.12 11.012 11.123
## 16745 1415672_at 11.432 11.417 11.442 11.555 11.561 11.318 11.461 11.561
## 17051 1415673_at 7.535 7.382 7.566 7.162 7.403 7.342 7.213 7.581
## 21715 1415674_a_at 9.757 9.972 9.269 9.873 9.354 9.918 9.459 9.655
## 11364 1415675_at 9.029 9.009 9.245 9.282 9.415 9.098 9.06 8.937
## BXD19 BXD2 BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29
## 18098 9.939 9.871 10.077 10.159 9.746 9.89 10.286 10.177 9.959
## 21104 10.922 10.802 10.988 10.969 11.104 10.979 10.905 11.013 11.071
## 16745 11.575 11.426 11.367 11.328 11.499 11.531 11.666 11.5 11.502
## 17051 7.551 7.368 7.251 7.373 7.408 7.473 7.238 7.424 7.475
## 21715 9.544 9.557 9.46 9.322 9.795 9.758 9.628 9.263 9.671
## 11364 8.995 9.104 9.123 9.086 9.087 9.04 8.844 9.173 8.951
## BXD31 BXD32 BXD33 BXD34 BXD36 BXD38 BXD39 BXD40 BXD42
## 18098 9.882 10.059 10.102 10.174 10.022 10.364 9.745 10.074 9.961
## 21104 10.996 10.854 11.084 11.059 10.923 11.053 11.03 11.067 10.891
## 16745 11.446 11.548 11.511 11.49 11.545 11.496 11.516 11.543 11.457
## 17051 7.491 7.419 7.386 7.496 7.369 7.265 7.13 7.149 7.383
## 21715 9.645 9.368 9.741 9.721 9.684 9.387 9.649 9.693 9.793
## 11364 9.131 9.108 9.122 9.153 9.008 9.138 9.126 9.233 9.06
## BXD5 BXD6 BXD8 BXD9
## 18098 10.16 10.069 9.956 10.142
## 21104 10.878 11.043 11.206 10.982
## 16745 11.374 11.504 11.531 11.7
## 17051 7.439 7.233 7.397 7.39
## 21715 9.84 9.519 9.34 9.16
## 11364 9.134 9.081 8.926 9.274
head(liver.gene.loc)
## probe_id Chromosome start_location end_location
## 1 1415670_at 6 87887971 87890759
## 3 1415671_at 8 105524469 105566040
## 4 1415672_at 8 23241325 23257080
## 5 1415673_at 5 129765557 129787253
## 6 1415674_a_at 9 44403758 44407548
## 7 1415675_at 2 32570857 32573571
dim(BXD.geno.loc)
## [1] 3811 3
dim(BXD.geno.SNP.eqtl)
## [1] 3811 31
head(BXD.geno.loc)
## Locus Chr position
## 1 rs6269442 1 3482275
## 2 rs6365999 1 4811062
## 3 rs6376963 1 5008089
## 4 rs3677817 1 5176058
## 5 rs8236463 1 5579193
## 6 rs6333200 1 6217921
head(BXD.geno.SNP.eqtl)
## Locus BXD1 BXD11 BXD12 BXD13 BXD14 BXD15 BXD16 BXD18 BXD19 BXD2
## 1 rs6269442 0 0 1 0 0 1 1 0 1 0
## 2 rs6365999 0 0 1 0 0 1 1 0 1 0
## 3 rs6376963 0 0 1 0 0 1 1 0 1 0
## 4 rs3677817 0 0 1 0 0 1 1 0 1 0
## 5 rs8236463 0 1 1 0 0 1 1 0 1 0
## 6 rs6333200 0 1 1 0 0 1 1 0 1 0
## BXD20 BXD21 BXD24 BXD24a BXD27 BXD28 BXD29 BXD31 BXD32 BXD33 BXD34 BXD36
## 1 1 1 0 0 0 1 1 0 0 0 0 0
## 2 1 1 0 0 0 1 1 0 0 0 0 0
## 3 1 1 0 0 0 1 1 0 0 0 0 0
## 4 1 1 0 0 0 1 1 0 0 0 0 0
## 5 1 1 0 0 0 1 1 0 0 0 0 0
## 6 1 1 0 0 0 1 1 0 0 0 0 0
## BXD38 BXD39 BXD40 BXD42 BXD5 BXD6 BXD8 BXD9
## 1 0 1 0 0 1 1 1 0
## 2 0 1 0 0 1 1 1 0
## 3 0 1 0 0 1 1 1 0
## 4 0 1 0 0 1 1 1 0
## 5 0 1 0 0 1 1 1 0
## 6 0 1 0 0 1 1 1 0