Nothing
allDuplicated <- function(x)
{
front <- duplicated(x)
back <- duplicated(x, fromLast = TRUE)
all_dup <- front | back
return(all_dup)
}
#' C version of g2c
#' @noRd
g2a.c <- function (g)
{
d <- 1 + 8 * (g - 1)
u <- 1 + ((1 + (sqrt(d) - 1) - 1) / 2)
u <- ceiling(u)
l = g - u * (u - 1) / 2
list (l=l,u=u)
}
is.miss <- function(data,data.int,miss.val=0)
{
if (data.int==2) # genotype
{
id <- array(FALSE, length(data))
for (i in 1:length(miss.val))
id <- id | data==miss.val[i]
}
else # allele
{
id <- array(FALSE, length(data[,1]))
for (i in 1:length(miss.val))
id <- id | apply(data==miss.val[i],1,any)
}
return (id)
}
#' Allele indices for a given haplotype ID in a multiallelic system
#' @section Usage:
#' revhap(loci,hapid)
#' @noRd
#'
#' @keywords internal
# adapted from gcp.c on 24/9/2004
revhap <- function(loci,hapid)
{
nloci <- length(loci)
nalleles <- vector("numeric",nloci)
nalleles[nloci] <- 1
m <- nloci
for(k in nloci:1)
{
m <- m - 1
nalleles[m] <- nalleles[m+1] * loci[k]
}
n <- length(hapid)
hap <- matrix(0,n,nloci)
for (i in 1:n)
{
l <- hapid[i] - 1
for (j in 1:nloci)
{
hap[i,j] <- floor(l/nalleles[j])
if (j==nloci) hap[i,j] <- l
else l <- l%%nalleles[j]
hap[i,j] <- hap[i,j] + 1
}
}
invisible(hap)
}
revhap.i <- function(loci,hapid)
{
nloci <- length(loci)
nalleles <- vector("numeric",nloci)
nalleles[nloci] <- 1
m <- nloci
for(k in nloci:1)
{
m <- m - 1
nalleles[m] <- nalleles[m+1] * loci[k]
}
hap <- vector("numeric",nloci)
l <- hapid - 1
for (j in 1:nloci)
{
hap[j] <- floor(l/nalleles[j])
if (j==nloci) hap[j] <- l
else l <- l%%nalleles[j]
hap[j] <- hap[j] + 1
}
invisible(hap)
}
#' as a2g
#' @noRd
gcode <- function(a1,a2) {
i <- ifelse(a1 < a2,a2,a1)
j <- ifelse(a1 < a2,a1,a2)
genocode <- i*(i-1)/2 + j
return(genocode)
}
#' Recovery of alleles from genotype(s)
#' @noRd
ungcode <- function(g) {
i <- 1 + floor((sqrt(8*g+1)-1)/2)
j <- g - i*(i-1)/2
# the following 2 lines were added as a patch to make this ungcode work:
i <- ifelse(j==0,i-1,i)
j <- ifelse(j==0,i,j)
return(cbind(j,i))
}
grec2g <- function (h, n, t)
{
hh <- h
for (i in 1:n)
{
hh[,i] <- t$alist[[i]]$allele[h[,i]]
}
invisible(hh)
}
#' an experimental function for PLEM format
#' @noRd
# a simple scheme to represent SNPs
# similar to a2g() and not unlike Mike Weale's Matlab function
# 10-2-2005 JHZ
m2plem <- function(a1,a2)
{
a <- vector("numeric")
# 1=A, 2=a
for (i in 1:length(a1))
{
if(a1[i]==1&a2[i]==2) {
a[i] <- 0
}
else if (a1[i]==1&a2[i]==1) {
a[i] <- 1
}
else if (a1[i]==2&a2[i]==2) {
a[i] <- 2
}
}
invisible(a)
}
#' PLEM format
#' @noRd
plem2m <- function(a)
{
a1 <- vector("numeric")
a2 <- vector("numeric")
# 1=A, 2=a
for (i in 1:length(a))
{
if(a[i]==0) {
a1[i] <- 1
a2[i] <- 2
}
else if (a[i]==1) {
a1[i] <- 1
a2[i] <- 1
}
else if (a[i]==2) {
a1[i] <- 2
a2[i] <- 2
}
}
invisible(list(a1,a2))
}
#' A function to combine imputation results
#' @noRd
micombine <- function (est, std.err, confidence = 0.95)
{
qstar <- est[[1]]
for (i in 2:length(est)) {
qstar <- cbind(qstar, est[[i]])
}
qbar <- apply(qstar, 1, mean)
u <- std.err[[1]]
for (i in 2:length(std.err)) {
u <- cbind(u, std.err[[i]])
}
if (!is.null(dimnames(qstar)[[1]]))
dimnames(u)[[1]] <- dimnames(qstar)[[1]]
u <- u^2
ubar <- apply(u, 1, mean)
bm <- apply(qstar, 1, var)
m <- dim(qstar)[2]
tm <- ubar + ((1 + (1/m)) * bm)
rem <- (1 + (1/m)) * bm/ubar
nu <- (m - 1) * (1 + (1/rem))^2
alpha <- 1 - (1 - confidence)/2
low <- qbar - qt(alpha, nu) * sqrt(tm)
up <- qbar + qt(alpha, nu) * sqrt(tm)
pval <- 2 * (1 - pt(abs(qbar/sqrt(tm)), nu))
fminf <- (rem + 2/(nu + 3))/(rem + 1)
result <- list(est = qbar, std.err = sqrt(tm), df = nu, signif = pval,
lower = low, upper = up, r = rem, fminf = fminf)
result
}
ReadGRMPLINK <- function(prefix,diag=1)
{
f <- paste(prefix,".genome",sep="")
g <- read.table(f,header=TRUE,as.is=TRUE)
L <- dim(g)[1]
N <- (1+sqrt(1+8*L))/2
pi_hat <- g["PI_HAT"]
PIHAT <- matrix(NA,N,N)
diag(PIHAT) <- diag
PIHAT[lower.tri(PIHAT)] <- pi_hat[,1]
PIHAT[upper.tri(PIHAT)] <- t(PIHAT)[upper.tri(PIHAT)]
idplink <- function (N)
{
id <- function(N, diag=TRUE)
{
irep <- function(i) rep(i,i)
iter <- function(i) 1:i
i <- unlist(lapply(1:N,irep))
j <- unlist(lapply(lapply(1:N,iter),cbind))
ij <- cbind(i,j)
if (diag) ij
else subset(ij,ij[,1]!=ij[,2])
}
ij <- id(N, diag=FALSE)
ij[order(ij[,2]),2:1]
}
fid1 <- g[1,"FID1"]
iid1 <- g[1,"IID1"]
fid <- c(fid1,g[1:(N-1),"FID2"])
iid <- c(iid1,g[1:(N-1),"IID2"])
ji <- idplink(N)
idn <- 1:N
ii <- cbind(i=idn,j=idn)
ij <- rbind(ji,ii)
idx <- ij[order(ij[,2]),2:1]
invisible(list(pihat=PIHAT[lower.tri(PIHAT,diag=TRUE)],PIHAT=PIHAT,id=cbind(fid,iid),idx=idx))
}
WriteGRMSAS <- function(grmlist,outfile="gwas")
{
for(p in c("foreign")) {
if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
if (!requireNamespace(p, quietly = TRUE))
warning(paste("WriteGRMSAS needs package `", p, "' to be fully functional; please install", sep=""))
}
}
with(grmlist,
{
parm <- 1
colnames(GRM) <- row <- paste("COL",1:dim(id)[1],sep="")
names(id) <- c("idn","id")
write.dta(data.frame(id,parm,row,GRM),paste(outfile,".dta",sep=""))
})
}
ReadGRMPCA <- function(prefix)
{
values <- scan(paste(prefix,".eigenval",sep=""))
vectors <- read.table(paste(prefix,".eigenvec",sep=""))
N <- dim(vectors)[1]
id <- vectors[,1:2]
vectors <- as.matrix(vectors[,-(1:2)])
s <- (values>0)
posGRM <- vectors[,s]%*%diag(values[s])%*%t(vectors[,s])
s <- (values<0)
negGRM <- vectors[,s]%*%diag(values[s])%*%t(vectors[,s])
invisible(list(N=N,posGRM=posGRM,negGRM=negGRM,id=id))
}
#' Variance of a ratio
#' @noRd
#'
#' @keywords internal
VR <- function(v1,vv1,v2,vv2,c12)
{
nv <- v2^2*vv1 + v1^2*vv2 - 2*v1*v2*c12
dv <- v2^4
nv/dv
}
lambda1000 <- function(lambda, ncases, ncontrols)
1 + (lambda - 1) * (1 / ncases + 1 / ncontrols)/( 1 / 1000 + 1 / 1000)
#' 3D Manhattan plot according to Sun, et al. (2018).
#'
#' @section Usage:
#' sun3d(xyz="INF1.merge.cis.vs.trans",
#' cols=c("id","chr1","pos1","chr2","pos2","gene","target","log10p","x","y","col"),
#' xy.scale=c(1.3e8,1.3e8),marker.size=4,log10p.max=400,
#' prefix=c("Sentinel","CHR","POS","CHR","POS","Gene","Target","-log10(p)"),
#' postfix="\u003c/br>",
#' json.file="d3.json",pretty=TRUE)
#' @noRd
#'
#' @keywords internal
sun3d <- function(xyz="INF1.merge.cis.vs.trans",
cols=c("id","chr1","pos1","chr2","pos2","gene","target","log10p","x","y","col"),
xy.scale=c(1.3e8,1.3e8),marker.size=4,log10p.max=400,
prefix=c("Sentinel","CHR","POS","CHR","POS","Gene","Target","-log10(p)"),
postfix="\u003c/br>",
json.file="d3.json",pretty=TRUE)
{
for(p in c("jsonlite", "plotly")) {
if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
if (!requireNamespace(p, quietly = TRUE))
warning(paste("sun3d needs package `", p, "' to be fully functional; please install", sep=""))
}
}
src <- list(
x = list(
layout = list(
margin = list(b = 40, l = 60, t = 25, r = 10),
scene = list(
xaxis = list(title = "pQTL position",
tickmode = "array",
autotick = FALSE,
tick0 = 1,
dtick = 1,
ticklen = 0,
tickwidth = 0,
tickfont = list(size = 10),
tickvals = as.list(1:24),
ticktext = as.list(c(1:22,"X","Y"))
),
yaxis = list(title = "Protein position",
tickmode = "array",
autotick = FALSE,
tick0 = 1,
dtick = 1,
ticklen = 0,
tickwidth = 0,
tickfont = list (size = 10),
tickvals = as.list(1:24),
ticktext = as.list(c(1:22,"X","Y"))
),
zaxis = list(title = "-log10(p)", tickfont = list(size = 10)),
camera = list(eye=list(x = -1.3, y = -1.2, z = 1.1)),
aspectmode = "manual",
aspectratio = list(x = 0.9, y = 1, z = 0.6)
),
legend = list(x = 10, y = 0.5),
xaxis = list(domain=list(0,1)),
yaxis = list(domain=list(0,1)),
title = "Scatterplot of sentinels",
showlegend = TRUE
),
source = "A",
config = list(
modeBarButtonsToAdd = list(name = "Collaborate"),
modeBarButtonsToRemove = list("sendDataToCloud")
),
data <- list(),
base_url = "https://plot.ly"
),
evals = list(),
jsHooks = list()
)
d <- read.csv(xyz,as.is=TRUE)
r <- qtl2dplot(d, plot=FALSE)
cuts <- with(r, abs(log10p) > log10p.max)
r <- within(r,{x=x/xy.scale[1]; y=y/xy.scale[2]; log10p[!cuts] <- abs(log10p[!cuts]); log10p[cuts] <- log10p.max})
fixes <- function(col,d) paste(paste(prefix[col],d[,col],sep=":"),postfix)
cistrans <- function(br,name,color)
{
t <- subset(r[cols],col==br)
cuts <- with(t, abs(log10p) == log10p.max)
n.cuts <- with(t,length(log10p[cuts]))
s <- rep('circle',nrow(t))
s[cuts] <- 'diamond'
list(x=t$x, y=t$y, z=t$log10p,
text=as.list(apply(sapply(c(1:8),fixes,t),1,paste,collapse=" ")),
type="scatter3d",
mode="markers",
name=name,
marker=list(color=color,symbol=s,size=marker.size))
}
cis <- cistrans("blue","cis","rgba(191,56,42,1)")
trans <- cistrans("red","trans","rgba(12,75,142,1)")
src$x$data <- list(cis,trans)
if(!is.null(json.file))
{
json <- jsonlite::toJSON(src,auto_unbox=TRUE,pretty=pretty)
sink(json.file)
print(json)
sink()
}
p <- plotly::plot_ly()
p <- with(src$x$data[[1]],plotly::add_trace(p, x=x, y=y, z=z, marker=marker, mode=mode, name=name, text=text, type=type))
p <- with(src$x$data[[2]],plotly::add_trace(p, x=x, y=y, z=z, marker=marker, mode=mode, name=name, text=text, type=type))
p <- with(src$x$layout, plotly::layout(p, scene=scene, xaxis=xaxis, yaxis=yaxis, margin=margin, title=title, showlegend=showlegend))
}
# https://plot.ly/r/reference/#scatter3d
# sed -i 's|<\\/br>|\\u003c/br>|g' d3.json
# plotly::toRGB( c('#BF382A', '#0C4B8E')) ==> "rgba(191,56,42,1)" "rgba(12,75,142,1)"
textbox <- function(label, name=NULL, gp=NULL, vp=NULL)
{
gt <- grid::gTree(label=label, name=name, gp=gp, vp=vp, cl="textboxtree")
grid::grid.draw(gt)
}
makeContent.textboxtree <- function(x)
{
t <- grid::textGrob(x$label, name="text")
rr <- grid::roundrectGrob(width=1.5*grid::grobWidth(t), height=1.5*grid::grobHeight(t), name="box")
grid::setChildren(x, grid::gList(t, rr))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.