Split a correlation matrix in blocks as independent as possible. This will find the splitting in blocks that minimize the sum of squared correlation between these blocks (i.e. everything outside these blocks).

snp_ldsplit(corr, thr_r2, min_size, max_size, max_K)

corr | Sparse correlation matrix. Usually, the output of |
---|---|

thr_r2 | Threshold under which squared correlations are ignored.
This is useful to avoid counting noise, which should give clearer patterns
of costs vs. number of blocks. It is therefore possible to have a splitting
cost of 0. If this parameter is used, then |

min_size | Minimum number of variants in each block. This is used not to have a disproportionate number of small blocks. |

max_size | Maximum number of variants in each block. This is used not to have blocks that are too large, e.g. to limit computational and memory requirements of applications that would use these blocks. For some long-range LD regions, it may be needed to allow for large blocks. |

max_K | Maximum number of blocks to consider. All optimal solutions for K
from 1 to |

A tibble with five columns:

`$n_block`

: Number of blocks.`$cost`

: The sum of squared correlations outside the blocks.`$perc_kept`

: Percentage of initial non-zero values kept within the blocks defined.`$block_num`

: Resulting block numbers for each variant.`$all_last`

: Last index of each block.`$all_size`

: Sizes of the blocks.

if (FALSE) { corr <- readRDS(url("https://www.dropbox.com/s/65u96jf7y32j2mj/spMat.rds?raw=1")) THR_R2 <- 0.01 (res <- snp_ldsplit(corr, thr_r2 = THR_R2, min_size = 10, max_size = 50, max_K = 50)) library(ggplot2) qplot(n_block, cost, data = res) + theme_bw(16) + scale_y_log10() all_ind <- head(res$all_last[[6]], -1) ## Transform sparse representation into (i,j,x) triplets corrT <- as(corr, "dgTMatrix") upper <- (corrT@i <= corrT@j & corrT@x^2 >= THR_R2) df <- data.frame( i = corrT@i[upper] + 1L, j = corrT@j[upper] + 1L, r2 = corrT@x[upper]^2 ) df$y <- (df$j - df$i) / 2 ggplot(df) + geom_point(aes(i + y, y, color = r2), size = rel(0.5)) + coord_fixed() + scale_color_gradientn(colours = rev(colorRamps::matlab.like2(100))) + theme_minimal() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) + geom_vline(xintercept = all_ind + 0.5, linetype = 3) + labs(x = "Position", y = NULL) + scale_alpha(guide = 'none') }