ld.snp takes an object of snp.matrix class and suitable
range and depth and calculation the pairwise D', $r^2$, LOD and return
the result as a snp.dprime object.
An object of snp.matrix class with M samples of
N snps
depth
The depth or lag of pair-wise calculation. Should be
between 1 and N-1; default 100. Using 0 (an invalid value) is the
same as picking the maximum
start
The index of the start of the range of interest. Should
be between 1 and (N-1); default 1
end
The index of the end of the range of interest. Should
be between 2 and N. default N.
signed.r
Boolean for whether to returned signed $r$ values
instead of $r^2$
Details
The cubic equation and quadratic equation solver code is borrowed from GSL (GNU Scientific Library).
Value
return a snp.dprime
object, which is a list of 3 named
matrices dprime, rsq2 (or r depending on the input), lod, and an attribute
snp.names for the list of snps involved. (Note that if $x$ snps
are involved, the row numbers of the 3 matrices are $(x-1)$).
Only one of rsq2 or r is present.
dprime
D'
rsq2
$r^2$
r
signed $r$
lod
Log of Odd's
All the matrices are defined such that the ($n, m$)th entry is the
pair-wise value between the ($n$)th snp and the $(n+m)$th snp. Hence
the lower right triangles are always filled with zeros. (See example
section for the actual layout)
Invalid values are represented by an out-of-range value - currently we
use -1 for D', $r^2$ (both of which are between 0 and 1), and -2 for
$r$ (valid values are between -1 and +1). lod is set to zero in
most of these invalid cases. (lod can be any value so it is not indicative).
Note
The output snp.dprime object is suitable
for input to plot.snp.dprime for drawing.
The speed of “ld.snp” LD calculation, on a single-processor
opteron 2.2GHz box:
unsigned $r^2$, 13191 snps, depth 100 = 36.4 s (~ 1.3 mil pairs)
signed r , 13191 snps, depth 100 = 40.94s (~ 1.3 mil pairs)
signed r , 13191 snps, depth 1500 = 582s (~ 18.5 mil pairs)
For depth=1500, it uses 500MB just for the three matrices. So I
actually cannot do the full depth at ~13,000; full depth should be
under 50 minutes for 87 mil pairs, even in the signed-r version.
The LD code can be ran outside of R - mainly for debugging:
Clayton, D.G. and Leung, Hin-Tak (2007) An R package for analysis of
whole-genome association studies.
Human Heredity64:45-51.
GSL (GNU Scientific Library) http://www.gnu.org/software/gsl/
See Also
snp.dprime-class,
plot.snp.dprime, ld.with
Examples
# LD stats between 500 SNPs at a depth of 50
data(testdata)
ldinfo <- ld.snp(Autosomes, start=1, end=500, depth=50)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(chopsticks)
Loading required package: survival
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/chopsticks/ld.snp.Rd_%03d_medium.png", width=480, height=480)
> ### Name: ld.snp
> ### Title: Function to calculate pairwise D', $r^2$
> ### Aliases: ld.snp
> ### Keywords: dplot htest models
>
> ### ** Examples
>
> # LD stats between 500 SNPs at a depth of 50
> data(testdata)
> ldinfo <- ld.snp(Autosomes, start=1, end=500, depth=50)
Information: The input contains 400 samples with 9445 snps
... Done
>
>
>
>
>
> dev.off()
null device
1
>