5.1 Pedigree drawing

I have included the original file for theR News as well as put examples in separatevignettes. They can be accessed viavignette("rnews",package="gap.examples") andvignette("pedtodot", package="gap.examples"), respectively.

We also demonstrate through pedigree 10081 example⁵ withpedtodot_verbatim.

library(gap)#> Loading required package: gap.datasets#> gap version 1.6knitr::kable(p3,caption="An example pedigree")

library(DOT)# one can see the diagram in RStudiopedtodot_verbatim(p3,run=TRUE,toDOT=TRUE,return="verbatim")#> Pedigree  10081 #> running DOT::dot on 10081library(DiagrammeR)pedtodot_verbatim(p3)#> Pedigree  10081grViz(readr::read_file("10081.dot"))

5.2 Kinship calculation

Next, I will provide an example for kinship calculation usingkin.morgan.It is recommended that individuals in a pedigree are ordered so that parentsalways precede their children. In this regard, packagepedigree can beused, and packagekinship2 can be used to produce pedigree diagram as withkinship matrix.

The pedigree diagram is as follows,

# pedigree diagramdata(lukas, package="gap.datasets")library(kinship2)#> Loading required package: Matrix#> Loading required package: quadprogped <- with(lukas,pedigree(id,father,mother,sex))plot(ped,cex=0.4)

Figure 5.2: A pedigree diagram

We then turn to the kinship calculation.

# unordered individualsgk1 <- kin.morgan(lukas)write.table(gk1$kin.matrix,"gap_1.txt",quote=FALSE)library(kinship2)kk1 <- kinship(lukas[,1],lukas[,2],lukas[,3])write.table(kk1,"kinship_1.txt",quote=FALSE)d <- gk1$kin.matrix-kk1sum(abs(d))#> [1] 2.443634# order individuals so that parents precede their childrenlibrary(pedigree)op <- orderPed(lukas)olukas <- lukas[order(op),]gk2 <- kin.morgan(olukas)write.table(olukas,"olukas.csv",quote=FALSE)write.table(gk2$kin.matrix,"gap_2.txt",quote=FALSE)kk2 <- kinship(olukas[,1],olukas[,2],olukas[,3])write.table(kk2,"kinship_2.txt",quote=FALSE)z <- gk2$kin.matrix-kk2sum(abs(z))#> [1] 0

We see that in the second case, the result agrees withkinship2.

6.1 Family-based design

The example is as follows,

options(width=150)library(gap)models <- matrix(c(         4.0, 0.01,         4.0, 0.10,         4.0, 0.50,          4.0, 0.80,         2.0, 0.01,         2.0, 0.10,         2.0, 0.50,         2.0, 0.80,         1.5, 0.01,             1.5, 0.10,         1.5, 0.50,         1.5, 0.80), ncol=2, byrow=TRUE)outfile <- "fbsize.txt"cat("gamma","p","Y","N_asp","P_A","H1","N_tdt","H2","N_asp/tdt",    "L_o","L_s\n",file=outfile,sep="\t")for(i in 1:12) {    g <- models[i,1]    p <- models[i,2]    z <- fbsize(g,p)    cat(z$gamma,z$p,z$y,z$n1,z$pA,z$h1,z$n2,z$h2,z$n3,        z$lambdao,z$lambdas,file=outfile,append=TRUE,sep="\t")    cat("\n",file=outfile,append=TRUE)}table1 <- read.table(outfile,header=TRUE,sep="\t")nc <- c(4,7,9)table1[,nc] <- ceiling(table1[,nc])dc <- c(3,5,6,8,10,11)table1[,dc] <- round(table1[,dc],2)unlink(outfile)knitr::kable(table1,caption="Power/Sample size of family-based designs")

As for APOE4 and Alzheimer’s¹⁴

g <- 4.5p <- 0.15alz <- data.frame(fbsize(g,p))knitr::kable(alz,caption="Power/Sample size of study on Alzheimer's disease")

6.2 Population-based design

The example is as follows,

library(gap)kp <- c(0.01,0.05,0.10,0.2)models <- matrix(c(          4.0, 0.01,          4.0, 0.10,          4.0, 0.50,           4.0, 0.80,          2.0, 0.01,          2.0, 0.10,          2.0, 0.50,          2.0, 0.80,          1.5, 0.01,              1.5, 0.10,          1.5, 0.50,          1.5, 0.80), ncol=2, byrow=TRUE)outfile <- "pbsize.txt"cat("gamma","p","p1","p5","p10","p20\n",sep="\t",file=outfile)for(i in 1:dim(models)[1]){   g <- models[i,1]   p <- models[i,2]   n <- vector()   for(k in kp) n <- c(n,ceiling(pbsize(k,g,p)))   cat(models[i,1:2],n,sep="\t",file=outfile,append=TRUE)   cat("\n",file=outfile,append=TRUE)} table5 <- read.table(outfile,header=TRUE,sep="\t")knitr::kable(table5,caption="Sample size of population-based design")

6.3 Case-cohort design

We obtain results for ARIC and EPIC studies.

library(gap)# ARIC studyoutfile <- "aric.txt"n <- 15792pD <- 0.03p1 <- 0.25alpha <- 0.05theta <- c(1.35,1.40,1.45)beta <- 0.2s_nb <- c(1463,722,468)cat("n","pD","p1","hr","q","power","ssize\n",file=outfile,sep="\t")for(i in 1:3){  q <- s_nb[i]/n  power <- ccsize(n,q,pD,p1,log(theta[i]),alpha,beta,TRUE)  ssize <- ccsize(n,q,pD,p1,log(theta[i]),alpha,beta,FALSE)  cat(n,"\t",pD,"\t",p1,"\t",theta[i],"\t",q,"\t",      signif(power,3),"\t",ssize,"\n",      file=outfile,append=TRUE)}read.table(outfile,header=TRUE,sep="\t")#>       n   pD   p1   hr          q power ssize#> 1 15792 0.03 0.25 1.35 0.09264184   0.8  1463#> 2 15792 0.03 0.25 1.40 0.04571935   0.8   722#> 3 15792 0.03 0.25 1.45 0.02963526   0.8   468unlink(outfile)# EPIC studyoutfile <- "epic.txt"n <- 25000alpha <- 0.00000005beta <- 0.2s_pD <- c(0.3,0.2,0.1,0.05)s_p1 <- seq(0.1,0.5,by=0.1)s_hr <- seq(1.1,1.4,by=0.1)cat("n","pD","p1","hr","alpha","ssize\n",file=outfile,sep="\t")# direct calculationfor(pD in s_pD){   for(p1 in s_p1)   {      for(hr in s_hr)      {         ssize <- ccsize(n,q,pD,p1,log(hr),alpha,beta,FALSE)         if (ssize>0) cat(n,"\t",pD,"\t",p1,"\t",hr,"\t",alpha,"\t",                          ssize,"\n",                          file=outfile,append=TRUE)      }   }}knitr::kable(read.table(outfile,header=TRUE,sep="\t"),caption="Sample size of case-cohort designs")

unlink(outfile)

7.1 Genome-wide association

The following code is used to obtain a Q-Q plot viaqqunif function,

library(gap)u_obs <- runif(1000)r <- qqunif(u_obs,pch=21,bg="blue",bty="n")u_exp <- r$yhits <- u_exp >= 2.30103points(r$x[hits],u_exp[hits],pch=21,bg="green")legend("topleft",sprintf("GC.lambda=%.4f",gc.lambda(u_obs)))

Figure 7.1: A Q-Q plot

Based on a chicken genome scan data, the code below generates a Manhattan plot, demonstratingthe use of gaps to separate chromosomes.

ord <- with(w4,order(chr,pos))w4 <- w4[ord,]oldpar <- par()par(cex=0.6)colors <- c(rep(c("blue","red"),15),"red")mhtplot(w4,control=mht.control(colors=colors,gap=1000),pch=19,srt=0)axis(2,cex.axis=2)suggestiveline <- -log10(3.60036E-05)genomewideline <- -log10(1.8E-06)abline(h=suggestiveline, col="blue")abline(h=genomewideline, col="green")abline(h=0)

Figure 7.2: A genome-wide association study on chickens

The code below obtains a Manhattan plot with gene annotation¹⁵,

data <- with(mhtdata,cbind(chr,pos,p))glist <- c("IRS1","SPRY2","FTO","GRIK3","SNED1","HTR1A","MARCH3","WISP3",           "PPP1R3B","RP1L1","FDFT1","SLC39A14","GFRA1","MC4R")hdata <- subset(mhtdata,gene%in%glist)[c("chr","pos","p","gene")]color <- rep(c("lightgray","gray"),11)glen <- length(glist)hcolor <- rep("red",glen)  par(las=2, xpd=TRUE, cex.axis=1.8, cex=0.4)ops <- mht.control(colors=color,yline=1.5,xline=3)hops <- hmht.control(data=hdata,colors=hcolor)mhtplot(data,ops,hops,pch=19)axis(2,pos=2,at=1:16,cex.axis=0.5)

Figure 7.3: A Manhattan plot with gene annotation

All these look familiar, so revised form of the function calledmhtplot2 wascreated for additional features such as centering the chromosome ticks, allowing formore sophisticated coloring schemes, using prespecified fonts, etc. Please refer tothe function’s documentation example.

We could also go further with a circos Manhattan plot,

circos.mhtplot(mhtdata, glist)#> Note: 11 points are out of plotting region in sector 'chr16', track '3'.

Figure 7.4: A circos Manhattan plot

and a version with y-axis,

require(gap.datasets)library(dplyr)#> #> Attaching package: 'dplyr'#> The following objects are masked from 'package:stats':#> #>     filter, lag#> The following objects are masked from 'package:base':#> #>     intersect, setdiff, setequal, uniontestdat <- mhtdata[c("chr","pos","p","gene","start","end")] %>%           rename(log10p=p) %>%           mutate(chr=paste0("chr",chr),log10p=-log10(log10p))dat <- mutate(testdat,start=pos,end=pos) %>%       select(chr,start,end,log10p)labs <- subset(testdat,gene %in% glist) %>%        group_by(gene,chr,start,end) %>%        summarize() %>%        mutate(cols="blue") %>%        select(chr,start,end,gene,cols)#> `summarise()` has grouped output by 'gene', 'chr', 'start'. You can override using the `.groups` argument.labs[2,"cols"] <- "red"ticks <- 0:2*5circos.mhtplot2(dat,labs,ticks=ticks,ymax=max(ticks))#> Note: 43 points are out of plotting region in sector 'chr16', track '5'.

Figure 7.5: Another circos Manhattan plot

As a side note, the data is used bymanhattanly.

#{r mhttest, fig.cap="Manhattanly plot", fig.height=8, fig.width=8, plotly=TRUE}library(manhattanly)mhttest <- manhattanly(mhtdata, chr = "chr", bp = "pos",                       snp = "rsn", annotation1 = "gene", suggestiveline = TRUE,                       annotation2 = "rsn", p = "p")mhttesthtmlwidgets::saveWidget(mhttest,"mhttest.html")

We now experiment with Miami plot,

mhtdata <- within(mhtdata,{pr=p})miamiplot(mhtdata,chr="chr",bp="pos",p="p",pr="pr",snp="rsn")

Figure 7.6: Miami plots

# An alternative implementationgwas <- select(mhtdata,chr,pos,p) %>%        mutate(z=qnorm(p/2))chrmaxpos <- miamiplot2(gwas,gwas,name1="Batch 2",name2="Batch 1",z1="z",z2="z")#> Warning in max(c(dat1$pos[dat1$chr == i], dat2$pos[dat2$chr == i]), na.rm = TRUE): no non-missing arguments to max; returning -InflabelManhattan(chr=c(2,16),pos=c(226814165,52373776),name=c("AnonymousGene","FTO"),gwas,gwasZLab="z",chrmaxpos=chrmaxpos)

Figure 7.7: Miami plots

and a truncated Manhattan plot, noting that only data points with -log10(P)>=2 are shown,

par(oma=c(0,0,0,0), mar=c(5,6.5,1,1))mhtplot.trunc(filter(IL.12B,log10P>=2), chr="Chromosome", bp="Position", z="Z",   snp="MarkerName",   suggestiveline=FALSE, genomewideline=-log10(5e-10),   cex.mtext=1.2, cex.text=1.2,   annotatelog10P=-log10(5e-10), annotateTop = FALSE,   highlight=with(genes,gene),   mtext.line=3, y.brk1=115, y.brk2=300, trunc.yaxis=TRUE, delta=0.01,   cex.axis=1.5, cex=0.8, font=3, font.axis=1.5,   y.ax.space=20,   col = c("blue4", "skyblue"))

Figure 7.8: Association of IL-12B

The code below obtains a regional association plot with theasplot function,

asplot(CDKNlocus, CDKNmap, CDKNgenes, best.pval=5.4e-8, sf=c(3,6))#> - CDKN2A #> - CDKN2Btitle("CDKN2A/CDKN2B Region")

Figure 7.9: A regional association plot

The function predates the currently popularlocuszoom software but leavesthe option open for generating such plots on the fly and locally.

Note that all these can serve as templates to customize features of your own.

7.2 Copy number variation

A plot of copy number variation (CNV) is shown here,

cnvplot(gap.datasets::cnv)

Figure 7.10: A CNV plot

7.3 Effect size plot

The code below obtains an effect size plot via theESplot function.

library(gap)rs12075 <- data.frame(id=c("CCL2","CCL7","CCL8","CCL11","CCL13","CXCL6","Monocytes"),                      b=c(0.1694,-0.0899,-0.0973,0.0749,0.189,0.0816,0.0338387),                      se=c(0.0113,0.013,0.0116,0.0114,0.0114,0.0115,0.00713386))ESplot(rs12075)

Figure 7.11: rs12075 and traits

7.4 Forest plot

It draws many forest plots given a list of variants, e.g.,

data(OPG,package="gap.datasets")meta::settings.meta(method.tau="DL")METAL_forestplot(OPGtbl,OPGall,OPGrsid,width=6.75,height=5,digits.TE=2,digits.se=2,                 col.diamond="black",col.inside="black",col.square="black")#> Joining with `by = join_by(MarkerName)`#> Joining with `by = join_by(MarkerName)`

Figure 7.12: Forest plots

Figure 7.13: Forest plots

METAL_forestplot(OPGtbl,OPGall,OPGrsid,package="metafor",method="FE",xlab="Effect",                 showweights=TRUE)#> Joining with `by = join_by(MarkerName)`#> Joining with `by = join_by(MarkerName)`

Figure 7.14: Forest plots

Figure 7.15: Forest plots

Figure 8.1: Normal(0,1) distribution

8.1 log(p) and log10(p)

First thing first, here are the answers for log(p) and log10(p) given z,

# log(p) for a standard normal deviate z based on log()logp <- function(z) log(2)+pnorm(-abs(z), lower.tail=TRUE, log.p=TRUE)# log10(p) for a standard normal deviate z based on log()log10p <- function(z) log(2, base=10)+pnorm(-abs(z), lower.tail=TRUE, log.p=TRUE)/log(10)

Notelogp() will be used for functions such asqnorm() as in functioncs() whereaslog10p() is more appropriate for Manhattan plot and used insentinels().

8.2 Rationale

We start with z=1.96 whose corresponding p value is approximately 0.05.

2*pnorm(-1.96,lower.tail=TRUE)

giving an acceptable value 0.04999579, so we proceed to get log10(p)

log10(2)+log10(pnorm(-abs(z),lower.tail=TRUE))

leading to the expression above from the fact that log10(X)=log(X)/log(10) since log(),being the natural log function, ln() – so log(exp(1)) = 1, in R, works far better onthe numerator of the second term. The use of -abs() just makes sure we are working onthe lower tail of the standard Normal distribution from which our p value is calculated.

8.3 Benchmark

Now we have a stress test,

z <- 20000-log10p(z)

giving -log10(p) = 86858901.

8.4 Multiple precision arithmetic

We would be curious about the p value itself as well, which is furnished with the Rmpfr package

require(Rmpfr)2*pnorm(mpfr(-abs(z),100),lower.tail=TRUE,log.p=FALSE)mpfr(log(2),100) + pnorm(mpfr(-abs(z),100),lower.tail=TRUE,log.p=TRUE)

giving p = 1.660579603192917090365313727164e-86858901 and -log(p) = -200000010.1292789076808554854177,respectively. To carry on we have -log10(p) = -log(p)/log(10)=86858901.

To make -log10(p) usable in R we obtain it directly through

as.numeric(-log10(2*pnorm(mpfr(-abs(z),100),lower.tail=TRUE)))

which actually yields exactly the same 86858901.

If we go very far to have z=50,000. then -log10(p)=542868107 but we have less luck with Rmpfr.

One may wonder the P value in this case, which is 6.6666145952e-542868108 or simply 6.67e-542868108.

The magic function for doing this is defined as follows,

pvalue <- function (z, decimals = 2){    lp <- -log10p(z)    exponent <- ceiling(lp)    base <- 10^-(lp - exponent)    paste0(round(base, decimals), "e", -exponent)}

and it is more appropriate to express p values in scientific format so they can behandled as follows,

log10pvalue <- function(p=NULL,base=NULL,exponent=NULL){  if(!is.null(p))  {    p <- format(p,scientific=TRUE)    p2 <- strsplit(p,"e")    base <- as.numeric(lapply(p2,"[",1))    exponent <- as.numeric(lapply(p2,"[",2))  } else if(is.null(base) | is.null(exponent)) stop("base and exponent should both be specified")  log10(base)+exponent}

used aslog10pvalue(p) when p<=1e-323, or log10pvalue(base=1,exponent=-323) otherwise.

One can also derive logpvalue for natural base (e) similarly.

We end with a quick look-up table

require(gap)zlist <- c(5,10,30,40,50,100,500,1000,2000,3000,5000)zp <- sapply(zlist,function(z) {c(z,pvalue(z),logp(z),log10p(z))})rownames(zp) <- c("z","P","log(P)","log10(P)")knitr::kable(t(zp),caption="z, P, log(P) and log10(P)")

8.5 Application

Themhtplot.trunc() function accepts three types of arguments:

16. P values of association statistics, which could be very small.

log10p. log10(P).

26. normal statistics that could be very large.

In all three cases, a log10(P) counterpart is obtained internally and to accommodate extreme value,the y-axis allows for truncation leaving out a given range to highlight the largest.

See the IL-12B example above.

9.1 Some preparations

Let\(\mbox{x} = SNP\ dosage\). Note that\(\mbox{Var}(\mbox{x})=2f(1-f)\),\(f=MAF\) or\(1-MAF\) by symmetry.

Our linear regression model is\(\mbox{y}=a + b\mbox{x} + e\). We have\(\mbox{Var}(\mbox{y}) = b^2\mbox{Var}(\mbox{x}) + \mbox{Var}(e)\). Moreover,\(\mbox{Var}(b)=\mbox{Var}(e)(\mbox{x}'\mbox{x})^{-1}=\mbox{Var}(e)/S_\mbox{xx}\), we have\(\mbox{Var}(e) = \mbox{Var}(b)S_\mbox{xx} = N \mbox{Var}(b) \mbox{Var}(\mbox{x})\). Consequently, let\(z = {b}/{SE(b)}\), we have

\[\begin{eqnarray*}\mbox{Var}(\mbox{y}) &=& \mbox{Var}(\mbox{x})(b^2+N\mbox{Var}(b)) \hspace{100cm} \cr &=& \mbox{Var}(\mbox{x})\mbox{Var}(\mbox{b})(z^2+N) \cr &=& 2f(1-f)(z^2+N)\mbox{Var}(b)\end{eqnarray*}\]

Moreover, the mean and the variance of the multiple correlation coefficient or the coefficient of determination (\(R^2\)) are known¹⁷ to be\({1}/{(N-1)}\) and\({2(N-2)}/{\left[(N-1)^2(N+1)\right]}\), respectively.

We also need some established results of a ratio (R/S)¹, i.e., the mean

\[\begin{align}E(R/S) \approx \frac{\mu_R}{\mu_S}-\frac{\mbox{Cov}(R,S)}{\mu_S^2}+\frac{\sigma_S^2\mu_R}{\mu_S^3} \hspace{100cm}\tag{9.1}\end{align}\]

and more importantly the variance

\[\begin{align}\mbox{Var}(R/S) \approx \frac{\mu_R^2}{\mu_S^2} \left[ \frac{\sigma_R^2}{\mu_R^2} -2\frac{\mbox{Cov}(R,S)}{\mu_R\;\mu_S} +\frac{\sigma_S^2}{\mu_S^2} \right] \hspace{100cm}\tag{9.2}\end{align}\]

where\(\mu_R\),\(\mu_S\),\(\sigma_R^2\),\(\sigma_S^2\) are the means and the variances for R and S, respectively.

Finally, we need some facts about\(\chi_1^2\),\(\chi^2\) distribution of one degree of freedom. For\(z \sim N(0,1)\),\(z^2\sim \chi_1^2\), whose mean and variance are 1 and 2, respectively.

We now have the following results.

9.2 Proportion of variance explained

We have

\[\begin{align}\mbox{PVE}_{\mbox{linear regression}} & = \frac{\mbox{Var}(b\mbox{x})}{\mbox{Var}(\mbox{y})} \hspace{100cm} \\ & = \frac{\mbox{Var}(\mbox{x})b^2}{\mbox{Var}(\mbox{x})(b^2+N\mbox{Var}(b))} \\ & = \frac{\mbox{z}^2}{\mbox{z}^2+N}\tag{9.3}\end{align}\]

On the other hand, for a simple linear regression\(R^2\equiv r^2\) where\(r\) is the Pearson correlation coefficient, which is readily from the\(t\)-statistic of the regression slope, i.e.,\(r={t}/{\sqrt{t^2+N-2}}\). so assuming\(t \equiv \ z \sim \chi_1^2\)

\[\begin{align}\mbox{PVE}_{t-\mbox{statistic}} & = \frac{\chi^2}{\chi^2+N-2} \hspace{100cm}\tag{9.4}\end{align}\]

To obtain coherent estimates of the asymptotic means and variances of both forms we resort to variance of a ratio (R/S). All the required elements are listed in a table below.

then we have the means and the variances for PVE.

Finally, our approximation of PVE for a protein with\(T\) independent pQTLs from the meta-analysis

Therefore they differ from the asymptotic results¹⁷ by ratios of\((N-2)(N+1)/(N-1)^2\) and\((N-2)/(N+1)\) forlinear regression and\(t\)-statistic, respectively.

Figure 9.1: Comparisons of Var(\(R^2\)) estimates

As the sample size increases, the estimates are quite close nevertheless quite small while the ratios approach 1 quickly from opposite sides after\(N\approx 300\).

9.3 Effect size and standard error

When\(\mbox{Var}(\mbox{y})=1\), as in cis eQTLGen¹⁸ data, we have\(b\) and its standard error (se) as follows,

\[\begin{align}b & = z/d \hspace{100cm} \\se & = 1/d\tag{9.5}\end{align}\]

where\(d = \sqrt{2f(1-f)(z^2+N)}\).

Now three functions are in place.

9.4 get_b_se

A record of the eQTLGen data is shown below

        SNP    Pvalue SNPChr  SNPPos AssessedAllele OtherAllele Zscore rs1003563 2.308e-06     12 6424577              A           G 4.7245            Gene GeneSymbol GeneChr GenePos NrCohorts NrSamples         FDR ENSG00000111321       LTBR      12 6492472        34     23991 0.006278872 BonferroniP hg19_chr hg19_pos AlleleA AlleleB allA_total allAB_total           1       12  6424577       A       G       2574        8483 allB_total AlleleB_all       7859   0.6396966

from which we obtain the effect size and its standard error as follows,

get_b_se(0.6396966,23991,4.7245)#>               b          se#> [1,] 0.04490488 0.009504684

9.5 get_pve_se

This function obtains proportion of explained variation (PVE) from n, z; itsstandard error is based on variance of the ratio (correction=TRUE) or\(r^2\).

9.6 get_sdy

We continue with the eQTLGen example above,

get_sdy(0.6396966,23991,0.04490488,0.009504684)#> [1] 1

and indeed the eQTLGen data were standardized.

10.1 Identification

Following an earlier implementation calledsentinels, a distance-based signalidentification based on GWAS summary statistics is avaialable fromqtlFinderfunction; to avoid the overhead of data-loading it works on a preselected listof variants from a GWAS and in this case a METAL output file.

10.2 Fine-mapping

We considered a region of interest (which could be approximately independent variants,e.g.,\(r^2 \le 0.5\)) using expressions that rely on effect sizes and their standarderrors¹⁹. More specifically, let Bayes factor (BF) for each variant in themeta-analysis be defined as\(ln(BF) \propto 0.5 \beta^2/SE^2\), where\(\beta\) and\(SE\)are the effect size and standard error from the meta-analysis, respectively. Theposterior probability (PP) for being causal for a particular variant is obtained as\(BF_i/\sum_{i=1}^TBF_i\), where\(i=1,\ldots,T\) indexes all variants considered in theregion. We generated credible sets within a given region by ranking all variants byPPs in descending order and calculating the number of variants required to reach acumulative probability of such as 99%.

The functioncs obtains credible set.

# zcat METAL/4E.BP1-1.tbl.gz | \# awk 'NR==1 || ($1==4 && $2 >= 187158034 - 1e6 && $2 < 187158034 + 1e6)' >  4E.BP1.ztbl <- within(read.delim("4E.BP1.z"),{logp <- logp(Effect/StdErr)})z <- cs(tbl)l <- cs(tbl,log_p="logp")

Note in particular that the implementation intends to avoid the naive summation inscenarios such as proteogenomic analysis containing exceptionally large BFs.

12.1 Themr function

The functionmr was originally developed to rework on data generated from GSMR²⁰, although it could be any harmonised data. The following example is from analysis of a real data on LIF-R protein and CVD/FEV1.

The MR analysis is as follows,

Figure 12.1: Mendelian randomization

Figure 12.2: Mendelian randomization

Figure 12.3: Mendelian randomization

Figure 12.4: Mendelian randomization

Figure 12.5: Mendelian randomization

Figure 12.6: Mendelian randomization

This is close to Ligthart et al.²¹ as used one time at workplace which turns to overlap with TwoSampleMR²².

12.2 Contrast of effect sizes

It would be of interest to contrast their effect sizes in the analysis above as well,

mr_names <- names(mrdat)LIF.R <- cbind(mrdat[grepl("SNP|LIF.R",mr_names)],trait="LIF.R"); names(LIF.R) <- c("SNP","b","se","trait")FEV1 <- cbind(mrdat[grepl("SNP|FEV1",mr_names)],trait="FEV1"); names(FEV1) <- c("SNP","b","se","trait")CAD <- cbind(mrdat[grepl("SNP|CAD",mr_names)],trait="CAD"); names(CAD) <- c("SNP","b","se","trait")mrdat2 <- within(rbind(LIF.R,FEV1,CAD),{y=b})library(ggplot2)p <- ggplot2::ggplot(mrdat2,aes(y = SNP, x = y))+     ggplot2::theme_bw()+     ggplot2::geom_point()+     ggplot2::facet_wrap(~ trait, ncol=3, scales="free_x")+     ggplot2::geom_segment(aes(x = b-1.96*se, xend = b+1.96*se, yend = SNP))+     ggplot2::geom_vline(lty=2, ggplot2::aes(xintercept=0), colour = 'red')+     ggplot2::xlab("Effect size")+     ggplot2::ylab("")p#> Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).#> Warning: Removed 1 row containing missing values or values outside the scale range (`geom_segment()`).

Figure 12.7: Combined forest plots for LIF.R, FEV1 and CAD

12.3 Forest plots

We illustrate withmr_forestplot(),

mr_forestplot(tnfb, colgap.forest.left="0.05cm", fontsize=14, leftlabs=c("Outcome","b","SE"),              common=FALSE, random=FALSE, print.I2=FALSE, print.pval.Q=FALSE, print.tau2=FALSE,              spacing=1.6,digits.TE=2,digits.se=2,xlab="Effect size",type.study="square",col.inside="black",col.square="black")

Figure 12.8: Forest plots for MR results on TNFB

mr_forestplot(tnfb, colgap.forest.left="0.05cm", fontsize=14,              leftcols="studlab", leftlabs="Outcome", plotwidth="3inch", sm="OR", rightlabs="ci",              sortvar=tnfb[["Effect"]],              common=FALSE, random=FALSE, print.I2=FALSE, print.pval.Q=FALSE, print.tau2=FALSE,              backtransf=TRUE, spacing=1.6,type.study="square",col.inside="black",col.square="black")

Figure 12.9: Forest plots for MR results on TNFB (no summary statistics)

mr_forestplot(tnfb,colgap.forest.left="0.05cm", fontsize=14,              leftcols=c("studlab"), leftlabs=c("Outcome"),              plotwidth="3inch", sm="OR", sortvar=tnfb[["Effect"]],              rightcols=c("effect","ci","pval"), rightlabs=c("OR","95%CI","P"),              digits=3, digits.pval=2, scientific.pval=TRUE,              common=FALSE, random=FALSE, print.I2=FALSE, print.pval.Q=FALSE, print.tau2=FALSE,              addrow=TRUE, backtransf=TRUE, spacing=1.6,type.study="square",col.inside="black",col.square="black")

Figure 12.10: Forest plots for MR results on TNFB (with P values)

13.1 chr_pos_a1_a2 and inv_chr_pos_a1_a2

They are functions to handle SNPid.

require(gap)s <- chr_pos_a1_a2(1,c(123,321),letters[1:2],letters[2:1])s#> [1] "chr1:123_A_B" "chr1:321_A_B"inv_chr_pos_a1_a2(s)#>               chr pos a1 a2#> chr1:123_A_B chr1 123  A  B#> chr1:321_A_B chr1 321  A  Binv_chr_pos_a1_a2("chr1:123-A_B",seps=c(":","-","_"))#>               chr pos a1 a2#> chr1:123-A_B chr1 123  A  B

13.2 ci2ms

Here is the documentation example on rs3784099 and breast cancer²³.

example(ci2ms)#> #> ci2ms> # rs3784099 and breast cancer recurrence/mortality#> ci2ms> ms <- ci2ms("1.28-1.72")#> #> ci2ms> print(ms)#> $m#> [1] 0.3945922#> #> $s#> [1] 0.07537491#> #> $direction#> [1] "+"#> #> #> ci2ms> # Vector input#> ci2ms> ci2 <- c("1.28-1.72","1.25-1.64")#> #> ci2ms> ms2 <- ci2ms(ci2)#> #> ci2ms> print(ms2)#> $m#> [1] 0.3945922 0.3589199#> #> $s#> [1] 0.07537491 0.06927492#> #> $direction#> [1] "+" "+"

13.3 gc.lambda

The definition is as follows,

gc.lambda <- function(x, logscale=FALSE, z=FALSE) {  v <- x[!is.na(x)]  n <- length(v)  if (z) {     obs <- v^2     exp <- qchisq(log(1:n/n),1,lower.tail=FALSE,log.p=TRUE)  } else {    if (!logscale)    {      obs <- qchisq(v,1,lower.tail=FALSE)      exp <- qchisq(1:n/n,1,lower.tail=FALSE)    } else {      obs <- qchisq(-log(10)*v,1,lower.tail=FALSE,log.p=TRUE)      exp <- qchisq(log(1:n/n),1,lower.tail=FALSE,log.p=TRUE)    }  }  lambda <- median(obs)/median(exp)  return(lambda)}# A simplified version is as follows,# obs <- median(chisq)# exp <- qchisq(0.5, 1) # 0.4549364# lambda <- obs/exp# see also estlambda from GenABEL and qq.chisq from snpStats# A related functionlambda1000 <- function(lambda, ncases, ncontrols)  1 + (lambda - 1) * (1 / ncases + 1 / ncontrols)/( 1 / 1000 + 1 / 1000)

13.4 invnormal

The function is widely used in various consortium analyses and defined as follows,

invnormal <- function(x) qnorm((rank(x,na.last="keep")-0.5)/sum(!is.na(x)))

An example use on data from Poisson distribution is as follows,

set.seed(12345)Ni <- rpois(50, lambda = 4); table(factor(Ni, 0:max(Ni)))#> #>  0  1  2  3  4  5  6  7  8  9 #>  2  4  6 11  8 10  4  2  2  1y <- invnormal(Ni)sd(y)#> [1] 0.9755074mean(y)#> [1] 0.002650512Ni <- 1:50y <- invnormal(Ni)mean(y)#> [1] -1.810184e-17sd(y)#> [1] 0.9973999

13.5 revStrand

This functions obtains allele(s) on the opposite strand.

alleles <- c("a","c","G","t")revStrand(alleles)#> [1] "t" "g" "C" "a"

13.6 snptest_sample

This is a function to output sample file for SNPTEST.

d <- data.frame(ID_1=1,ID_2=1,missing=0,PC1=1,PC2=2,D1=1,P1=10)snptest_sample(d,C=paste0("PC",1:2),D=paste0("D",1:1),P=paste0("P",1:1))

The commands above generates a file namedsnptest.sample.