### Example: Calculate GRM and display Source: https://privefl.github.io/bigsnpr/reference/bed_tcrossprodSelf.html This example demonstrates how to load a bed file, compute the GRM using `bed_tcrossprodSelf`, and then display a portion of the resulting matrix normalized by the number of columns. ```R bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") obj.bed <- bed(bedfile) K <- bed_tcrossprodSelf(obj.bed) K[1:4, 1:6] / ncol(obj.bed) ``` -------------------------------- ### Get Example Bedfile Source: https://privefl.github.io/bigsnpr/articles/demo.html Retrieves the example BED file from the bigsnpr package. This is the starting point for demonstrating genotype data handling. ```R # Get the example bedfile from package bigsnpr bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") ``` -------------------------------- ### Attach example bigSNP files Source: https://privefl.github.io/bigsnpr/reference/snp_attachExtdata.html Use this function to attach example bigSNP files. By default, it attaches 'example.bed'. You can also specify 'example-missing.bed' to test handling of missing data. ```R snp_attachExtdata(bedfile = c("example.bed", "example-missing.bed")) ``` -------------------------------- ### Example: Creating and Reading PLINK Files Source: https://privefl.github.io/bigsnpr/reference/snp_writeBed.html This example demonstrates how to create a fake bigSNP object, write it to PLINK format using snp_writeBed, and then read it back using snp_readBed. It verifies that the content of the original and newly created PLINK files are identical. ```R N <- 17 M <- 911 fake <- snp_fake(N, M) G <- fake$genotypes G[] <- sample(as.raw(0:3), size = length(G), replace = TRUE) # Write the object as a bed/bim/fam object tmp <- tempfile(fileext = ".bed") bed <- snp_writeBed(fake, tmp) # Read this new file for the first time rds <- snp_readBed(bed, backingfile = tempfile()) # Attach object in R session fake2 <- snp_attach(rds) # Same content all.equal(fake$genotypes[], fake2$genotypes[]) #> [1] TRUE all.equal(fake$fam, fake2$fam) #> [1] TRUE all.equal(fake$map, fake2$map) #> [1] TRUE # Two different backingfiles fake$genotypes$backingfile #> [1] "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file643850567c45.bk" fake2$genotypes$backingfile #> [1] "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file6438357c5f9b.bk" ``` -------------------------------- ### Example: Reading and attaching a bigSNP object Source: https://privefl.github.io/bigsnpr/reference/snp_readBed.html This example demonstrates how to read a bedfile using `snp_readBed`, store the data in a temporary directory, and then attach the bigSNP object using `snp_attach` for further analysis. It also shows how to inspect the structure and a subset of the genotype data. ```R (bedfile <- system.file("extdata", "example.bed", package = "bigsnpr")) #> [1] "C:/Users/au639593/AppData/Local/Temp/Rtmp8isrlL/temp_libpath2fc45f76188/bigsnpr/extdata/example.bed" # Reading the bedfile and storing the data in temporary directory rds <- snp_readBed(bedfile, backingfile = tempfile()) # Loading the data from backing files test <- snp_attach(rds) str(test) #> List of 3 #> $ genotypes:Reference class 'FBM.code256' [package "bigstatsr"] with 16 fields #> ..$ extptr : #> ..$ extptr_rw : #> ..$ nrow : int 517 #> ..$ ncol : int 4542 #> ..$ type : Named int 1 #> .. ..- attr(*, "names")= chr "unsigned char" #> ..$ backingfile : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file64381fb568d5.bk" #> ..$ is_read_only: logi FALSE #> ..$ address : #> ..$ address_rw : #> ..$ bk : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file64381fb568d5.bk" #> ..$ rds : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file64381fb568d5.rds" #> ..$ is_saved : logi TRUE #> ..$ type_chr : chr "unsigned char" #> ..$ type_size : int 1 #> ..$ file_size : num 2348214 #> ..$ code256 : num [1:256] 0 1 2 NA NA NA NA NA NA NA ... #> ..and 26 methods, of which 12 are possibly relevant: #> .. add_columns, as.FBM, bm, bm.desc, check_dimensions, #> .. check_write_permissions, copy#envRefClass, initialize, initialize#FBM, #> .. save, show#FBM, show#envRefClass #> $ fam :'data.frame': 517 obs. of 6 variables: #> ..$ family.ID : chr [1:517] "POP1" "POP1" "POP1" "POP1" ... #> ..$ sample.ID : chr [1:517] "IND0" "IND1" "IND2" "IND3" ... #> ..$ paternal.ID: int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ maternal.ID: int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ sex : int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ affection : int [1:517] 1 1 2 1 1 1 1 1 1 1 ... #> $ map :'data.frame': 4542 obs. of 6 variables: #> ..$ chromosome : int [1:4542] 1 1 1 1 1 1 1 1 1 1 ... #> ..$ marker.ID : chr [1:4542] "SNP0" "SNP1" "SNP2" "SNP3" ... #> ..$ genetic.dist: int [1:4542] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ physical.pos: int [1:4542] 112 1098 2089 3107 4091 5091 6107 7103 8090 9074 ... #> ..$ allele1 : chr [1:4542] "A" "T" "T" "T" ... #> ..$ allele2 : chr [1:4542] "T" "A" "A" "A" ... #> - attr(*, "class")= chr "bigSNP" dim(G <- test$genotypes) #> [1] 517 4542 G[1:8, 1:8] #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] #> [1,] 0 0 2 0 1 1 0 2 #> [2,] 1 0 1 0 0 1 0 2 #> [3,] 0 1 1 0 2 1 0 2 #> [4,] 0 0 2 0 2 1 0 2 #> [5,] 1 0 0 0 2 2 1 0 #> [6,] 0 1 0 0 2 1 0 0 #> [7,] 0 1 1 0 2 2 1 1 #> [8,] 1 0 1 1 1 1 0 1 ``` -------------------------------- ### Example: Get a table of related pairs without creating new bed files Source: https://privefl.github.io/bigsnpr/reference/snp_plinkKINGQC.html This example shows how to obtain a table of related pairs by setting `make.bed = FALSE`. The function returns a data frame containing kinship coefficients for related pairs, which can then be further analyzed. ```R df_rel <- snp_plinkKINGQC(plink2, bedfile, make.bed = FALSE, ncores = 2) str(df_rel) ``` -------------------------------- ### Example of attaching and inspecting a bigSNP Source: https://privefl.github.io/bigsnpr/reference/snp_subset.html This example demonstrates how to attach an external `bigSNP` file and inspect its structure using `str()`. It shows the components of the `bigSNP` object, including genotypes, family information, and marker information. ```r str(test <- snp_attachExtdata()) ``` -------------------------------- ### Example: Loading and inspecting a bigSNP object Source: https://privefl.github.io/bigsnpr/reference/snp_attach.html This example demonstrates the complete workflow of reading a bedfile, saving it to temporary backing files using `snp_readBed`, and then attaching it back into R using `snp_attach`. It concludes by showing the structure and dimensions of the loaded bigSNP object and its genotype matrix. ```r (bedfile <- system.file("extdata", "example.bed", package = "bigsnpr")) #> [1] "C:/Users/au639593/AppData/Local/Temp/Rtmp8isrlL/temp_libpath2fc45f76188/bigsnpr/extdata/example.bed" # Reading the bedfile and storing the data in temporary directory rds <- snp_readBed(bedfile, backingfile = tempfile()) # Loading the data from backing files test <- snp_attach(rds) str(test) #> List of 3 #> $ genotypes:Reference class 'FBM.code256' [package "bigstatsr"] with 16 fields #> ..$ extptr : #> ..$ extptr_rw : #> ..$ nrow : int 517 #> ..$ ncol : int 4542 #> ..$ type : Named int 1 #> .. ..- attr(*, "names")= chr "unsigned char" #> ..$ backingfile : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file643867a06069.bk" #> ..$ is_read_only: logi FALSE #> ..$ address : #> ..$ address_rw : #> ..$ bk : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file643867a06069.bk" #> ..$ rds : chr "C:\\Users\\au639593\\AppData\\Local\\Temp\\RtmpsvW07Y\\file643867a06069.rds" #> ..$ is_saved : logi TRUE #> ..$ type_chr : chr "unsigned char" #> ..$ type_size : int 1 #> ..$ file_size : num 2348214 #> ..$ code256 : num [1:256] 0 1 2 NA NA NA NA NA NA NA ... #> ..and 26 methods, of which 12 are possibly relevant: #> .. add_columns, as.FBM, bm, bm.desc, check_dimensions, #> .. check_write_permissions, copy#envRefClass, initialize, initialize#FBM, #> .. save, show#FBM, show#envRefClass #> $ fam :'data.frame': 517 obs. of 6 variables: #> ..$ family.ID : chr [1:517] "POP1" "POP1" "POP1" "POP1" ... #> ..$ sample.ID : chr [1:517] "IND0" "IND1" "IND2" "IND3" ... #> ..$ paternal.ID: int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ maternal.ID: int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ sex : int [1:517] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ affection : int [1:517] 1 1 2 1 1 1 1 1 1 1 ... #> $ map :'data.frame': 4542 obs. of 6 variables: #> ..$ chromosome : int [1:4542] 1 1 1 1 1 1 1 1 1 1 ... #> ..$ marker.ID : chr [1:4542] "SNP0" "SNP1" "SNP2" "SNP3" ... #> ..$ genetic.dist: int [1:4542] 0 0 0 0 0 0 0 0 0 0 ... #> ..$ physical.pos: int [1:4542] 112 1098 2089 3107 4091 5091 6107 7103 8090 9074 ... #> ..$ allele1 : chr [1:4542] "A" "T" "T" "T" ... #> ..$ allele2 : chr [1:4542] "T" "A" "A" "A" ... #> - attr(*, "class")= chr "bigSNP" dim(G <- test$genotypes) #> [1] 517 4542 G[1:8, 1:8] #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] #> [1,] 0 0 2 0 1 1 0 2 #> [2,] 1 0 1 0 0 1 0 2 #> [3,] 0 1 1 0 2 1 0 2 #> [4,] 0 0 2 0 2 1 0 2 #> [5,] 1 0 0 0 2 2 1 0 #> [6,] 0 1 0 0 2 1 0 0 #> [7,] 0 1 1 0 2 2 1 1 #> [8,] 1 0 1 1 1 1 0 1 ``` -------------------------------- ### Example: Generate new bed files after pruning related individuals Source: https://privefl.github.io/bigsnpr/reference/snp_plinkKINGQC.html This example demonstrates how to use `snp_plinkKINGQC` to generate new bed files with related individuals pruned. It first downloads PLINK 2, then calls the function with specified output file and number of cores. ```R bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") plink2 <- download_plink2(AVX2 = FALSE) bedfile2 <- snp_plinkKINGQC(plink2, bedfile, bedfile.out = tempfile(fileext = ".bed"), ncores = 2) ``` -------------------------------- ### Install bigstatsr and bigsnpr from GitHub Source: https://privefl.github.io/bigsnpr-extdoc/introduction.html Install the latest development versions of the packages from GitHub using the remotes package. ```r # install.packages("remotes") remotes::install_github("privefl/bigstatsr") remotes::install_github("privefl/bigsnpr") ``` -------------------------------- ### Example: Constructing Fake Genotype Data and Computing MAX3 Source: https://privefl.github.io/bigsnpr/reference/snp_MAX3.html This example demonstrates how to construct a fake genotype big.matrix, specify case/control phenotypes, and then compute MAX3 statistics using `snp_MAX3`. It also shows how to inspect the output and use `snp_qq` for visualizing p-value distribution. ```r set.seed(1) # constructing a fake genotype big.matrix N <- 50; M <- 1200 fake <- snp_fake(N, M) G <- fake$genotypes G[] <- sample(as.raw(0:3), size = length(G), replace = TRUE) G[1:8, 1:10] #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] #> [1,] 0 0 NA 0 2 1 0 NA NA 2 #> [2,] NA 0 NA 2 NA 2 2 2 2 1 #> [3,] 2 2 2 2 1 2 1 0 2 0 #> [4,] NA 1 0 NA 0 2 2 0 0 2 #> [5,] 0 0 NA NA 2 NA 2 2 NA NA #> [6,] 1 0 NA 2 2 NA 1 2 2 0 #> [7,] 2 NA 1 1 1 NA 1 2 1 NA #> [8,] 2 2 1 0 1 NA NA 2 0 1 # Specify case/control phenotypes fake$fam$affection <- rep(1:2, each = N / 2) # Get MAX3 statistics y01 <- fake$fam$affection - 1 str(test <- snp_MAX3(fake$genotypes, y01.train = y01)) #> Classes 'mhtest' and 'data.frame': 1200 obs. of 1 variable: #> $ score: num 1.7124 0.0178 2.7611 0.1309 0.4396 ... #> - attr(*, "transfo")=function (x) #> - attr(*, "predict")=function (xtr) # p-values are not well calibrated snp_qq(test) # genomic control is not of much help snp_qq(snp_gc(test)) # Armitage trend test (well calibrated because only one test) test2 <- snp_MAX3(fake$genotypes, y01.train = y01, val = 0.5) snp_qq(test2) ``` -------------------------------- ### Download Tutorial Data with runonce Source: https://privefl.github.io/bigsnpr-extdoc/index.html Use the runonce::download_file function to download necessary data files for tutorials. Ensure you have R and RStudio installed. ```r runonce::download_file( "https://figshare.com/ndownloader/files/38019072", dir = "tmp-data", fname = "GWAS_data.zip") # 109 MB runonce::download_file( "https://figshare.com/ndownloader/files/38019027", dir = "tmp-data", fname = "ref_freqs.csv.gz") # 46 MB runonce::download_file( "https://figshare.com/ndownloader/files/38019024", dir = "tmp-data", fname = "projection.csv.gz") # 44 MB runonce::download_file( "https://figshare.com/ndownloader/files/38077323", dir = "tmp-data", fname = "sumstats_CAD_tuto.csv.gz") # 16 MB runonce::download_file( "https://figshare.com/ndownloader/files/38247288", dir = "tmp-data", fname = "gen_pos_tuto.rds") # 2.5 MB ``` -------------------------------- ### Install bigsnpr from GitHub Source: https://privefl.github.io/bigsnpr/index.html Installs the bigsnpr package directly from its GitHub repository using the remotes package. This is useful for obtaining the latest development version. ```R # install.packages("remotes") remotes::install_github("privefl/bigsnpr") ``` -------------------------------- ### Install bigstatsr and bigsnpr from CRAN Source: https://privefl.github.io/bigsnpr-extdoc/introduction.html Use install.packages() to install the stable versions of both packages from CRAN. ```r install.packages("bigstatsr") install.packages("bigsnpr") ``` -------------------------------- ### Attach a "bigSNP" for examples and tests Source: https://privefl.github.io/bigsnpr/reference/index.html Attaches a bigSNP object specifically for use in examples and tests. ```APIDOC ## snp_attachExtdata() ### Description Attaches a "bigSNP" for examples and tests. ### Method Not applicable (R function) ### Endpoint Not applicable (R function) ### Parameters #### Path Parameters None #### Query Parameters None #### Request Body None ### Request Example ```R snp_attachExtdata("example_data_name") ``` ### Response #### Success Response (200) Returns a bigSNP object for examples/tests. #### Response Example ```R # Example of a bigSNP object structure for examples ``` ``` -------------------------------- ### Download BGEN files for conversion Source: https://privefl.github.io/bigsnpr-extdoc/inputs-and-formats.html Downloads example BGEN, BGI, and sample files required for converting BGEN format data. ```R bgen <- runonce::download_file( "https://enkre.net/cgi-bin/code/bgen/raw/3ec770a829a753282b5cb45afc3f4eda036b246705b76f9037b6cc98c41a4194?at=example.8bits.bgen", fname = "example.bgen") bgi <- runonce::download_file( "https://enkre.net/cgi-bin/code/bgen/raw/dc7276e0f0e2e096f58d2dac645aa5711de2cd64c3b29a07a80575e175344f78?at=example.8bits.bgen.bgi", fname = "example.bgen.bgi") sample <- runonce::download_file( "https://enkre.net/cgi-bin/code/bgen/raw/a3c4d8e4c132048a502dc00a3e51362f98eda5a2889df695ba260dc48c327fd9?at=example.sample", fname = "example.sample") ``` -------------------------------- ### Example: PRS Calculation and Visualization Source: https://privefl.github.io/bigsnpr/reference/snp_PRS.html This example demonstrates a full workflow: attaching external data, performing PCA for covariables, training a logistic regression model for GWAS, clumping SNPs, calculating PRS with various thresholds, and visualizing AUC performance. ```R test <- snp_attachExtdata() G <- big_copy(test$genotypes, ind.col = 1:1000) CHR <- test$map$chromosome[1:1000] POS <- test$map$physical.position[1:1000] y01 <- test$fam$affection - 1 # PCA -> covariables obj.svd <- snp_autoSVD(G, infos.chr = CHR, infos.pos = POS) #> Discarding 0 variant with MAC < 10 or MAF < 0.02. #> #> Phase of clumping (on MAF) at r^2 > 0.2.. keep 952 variants. #> #> Iteration 1: #> Computing SVD.. #> 26 outlier variants detected.. #> #> Iteration 2: #> Computing SVD.. #> 0 outlier variant detected.. #> #> Converged! # train and test set ind.train <- sort(sample(nrow(G), 400)) ind.test <- setdiff(rows_along(G), ind.train) # 117 # GWAS gwas.train <- big_univLogReg(G, y01.train = y01[ind.train], ind.train = ind.train, covar.train = obj.svd$u[ind.train, ]) # clumping ind.keep <- snp_clumping(G, infos.chr = CHR, ind.row = ind.train, S = abs(gwas.train$score)) # -log10(p-values) and thresolding summary(lpS.keep <- -predict(gwas.train)[ind.keep]) #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 0.000569 0.146282 0.344033 0.474463 0.640888 6.191733 thrs <- seq(0, 4, by = 0.5) nb.pred <- sapply(thrs, function(thr) sum(lpS.keep > thr)) # PRS prs <- snp_PRS(G, betas.keep = gwas.train$estim[ind.keep], ind.test = ind.test, ind.keep = ind.keep, lpS.keep = lpS.keep, thr.list = thrs) # AUC as a function of the number of predictors aucs <- apply(prs, 2, AUC, target = y01[ind.test]) library(ggplot2) #> Warning: package 'ggplot2' was built under R version 4.2.3 qplot(nb.pred, aucs) + geom_line() + scale_x_log10(breaks = nb.pred) + labs(x = "Number of predictors", y = "AUC") + theme_bigstatsr() #> Warning: `qplot()` was deprecated in ggplot2 3.4.0. ``` -------------------------------- ### Example: Attaching External Data and Counting Families Source: https://privefl.github.io/bigsnpr/reference/snp_getSampleInfos.html This example demonstrates how to attach external data using snp_attachExtdata() and then count the occurrences of family IDs within the attached data. This is a common first step before using snp_getSampleInfos. ```R test <- snp_attachExtdata() table(test$fam$family.ID) ``` -------------------------------- ### Example: Performing and Visualizing Fast Imputation Source: https://privefl.github.io/bigsnpr/reference/snp_fastImpute.html This example demonstrates how to use snp_fastImpute to impute missing genotypes and then visualize the results. It includes steps for attaching external data, performing imputation, and plotting the proportion of missing values against estimated imputation errors. Note that the imputation is not permanent until the modified genotype data is saved. ```R fake <- snp_attachExtdata("example-missing.bed") G <- fake$genotypes CHR <- fake$map$chromosome infos <- snp_fastImpute(G, CHR) infos[, 1:5] # Still missing values big_counts(G, ind.col = 1:10) # You need to change the code of G # To make this permanent, you need to save (modify) the file on disk fake$genotypes$code256 <- CODE_IMPUTE_PRED fake <- snp_save(fake) big_counts(fake$genotypes, ind.col = 1:10) # Plot for post-checking ## Here there is no SNP with more than 1% error (estimated) pvals <- c(0.01, 0.005, 0.002, 0.001); colvals <- 2:5 df <- data.frame(pNA = infos[1, ], pError = infos[2, ]) # base R plot(subset(df, pNA > 0.001), pch = 20) idc <- lapply(seq_along(pvals), function(i) { curve(pvals[i] / x, from = 0, lwd = 2, col = colvals[i], add = TRUE) }) legend("topright", legend = pvals, title = "p(NA & Error)", col = colvals, lty = 1, lwd = 2) # ggplot2 library(ggplot2) Reduce(function(p, i) { p + stat_function(fun = function(x) pvals[i] / x, color = colvals[i]) }, x = seq_along(pvals), init = ggplot(df, aes(pNA, pError))) + geom_point() + coord_cartesian(ylim = range(df$pError, na.rm = TRUE)) + theme_bigstatsr() } ``` -------------------------------- ### Prepare data and run GWAS for Q-Q plot Source: https://privefl.github.io/bigsnpr/reference/snp_qq.html This example demonstrates setting up genotype data, generating phenotype data, and performing a linear regression GWAS using `big_univLinReg`. The results are then used for plotting. ```R set.seed(9) test <- snp_attachExtdata() G <- test$genotypes y <- rnorm(nrow(G)) gwas <- big_univLinReg(G, y) ``` -------------------------------- ### Get sample information Source: https://privefl.github.io/bigsnpr/reference/index.html Retrieves sample information from a bigSNP object. ```APIDOC ## snp_getSampleInfos() ### Description Retrieves sample information from a bigSNP object. ### Method Not applicable (R function) ### Endpoint Not applicable (R function) ### Parameters #### Path Parameters None #### Query Parameters None #### Request Body None ### Request Example ```R snp_getSampleInfos(bigSNP_object) ``` ### Response #### Success Response (200) Returns a data frame with sample information. #### Response Example ```R # Sample information data frame ``` ``` -------------------------------- ### Example Usage of bed_cprodVec Source: https://privefl.github.io/bigsnpr/reference/bed_cprodVec.html Demonstrates how to use the bed_cprodVec function by loading a bed file, creating a sample vector, and then computing the cross-product. Shows the structure of the output. ```R bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") obj.bed <- bed(bedfile) y.row <- rep(1, nrow(obj.bed)) str(bed_cprodVec(obj.bed, y.row)) #> num [1:4542] 354 213 245 191 472 368 132 497 498 481 ... ``` -------------------------------- ### Example usage of bed_randomSVD Source: https://privefl.github.io/bigsnpr/reference/bed_randomSVD.html Demonstrates how to load a .bed file and perform a randomized partial SVD using `bed_randomSVD`. The output structure of the SVD result is also shown. ```R bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") obj.bed <- bed(bedfile) str(bed_randomSVD(obj.bed)) ``` -------------------------------- ### Example: Estimating ancestry proportions from GWAS data Source: https://privefl.github.io/bigsnpr/reference/snp_ancestry_summary.html This example demonstrates how to estimate ancestry proportions for a GWAS dataset. It involves downloading and processing summary statistics, reference frequencies, and projection data, then using `snp_match` and `snp_ancestry_summary` to calculate the proportions. Finally, it groups and sums the ancestry coefficients for different geographical regions. ```R if (FALSE) { # GWAS summary statistics for Epilepsy (supposedly in EUR+EAS+AFR) gz <- runonce::download_file( "http://www.epigad.org/gwas_ilae2018_16loci/all_epilepsy_METAL.gz", dir = "tmp-data") readLines(gz, n = 3) library(dplyr) sumstats <- bigreadr::fread2( gz, select = c("CHR", "BP", "Allele2", "Allele1", "Freq1"), col.names = c("chr", "pos", "a0", "a1", "freq") ) %>% mutate_at(3:4, toupper) # It is a good idea to filter for similar per-variant N (when available..) all_freq <- bigreadr::fread2( runonce::download_file("https://figshare.com/ndownloader/files/31620968", dir = "tmp-data", fname = "ref_freqs.csv.gz")) projection <- bigreadr::fread2( runonce::download_file("https://figshare.com/ndownloader/files/31620953", dir = "tmp-data", fname = "projection.csv.gz")) matched <- snp_match( mutate(sumstats, chr = as.integer(chr), beta = 1), all_freq[1:5], return_flip_and_rev = TRUE ) %>% mutate(freq = ifelse(`_REV_`, 1 - freq, freq)) res <- snp_ancestry_summary( freq = matched$freq, info_freq_ref = all_freq[matched$`_NUM_ID_`, -(1:5)], projection = projection[matched$`_NUM_ID_`, -(1:5)], correction = c(1, 1, 1, 1.008, 1.021, 1.034, 1.052, 1.074, 1.099, 1.123, 1.15, 1.195, 1.256, 1.321, 1.382, 1.443) ) # Some ancestry groups are very close to each other, and should be merged group <- colnames(all_freq)[-(1:5)] group[group %in% c("Scandinavia", "United Kingdom", "Ireland")] <- "Europe (North West)" group[group %in% c("Europe (South East)", "Europe (North East)")] <- "Europe (East)" tapply(res, factor(group, unique(group)), sum) } ``` -------------------------------- ### Example: Adding and Saving a 'p-values' Slot Source: https://privefl.github.io/bigsnpr/reference/snp_save.html This example demonstrates how to add a new slot (e.g., 'p-values') to the map data of a `bigSNP` object and then save these modifications using `snp_save`. It shows the state of the map data before saving, after re-reading without saving, and after re-reading with saving. ```R set.seed(1) # Reading example test <- snp_attachExtdata() # I can add whatever I want to an S3 class test$map$`p-values` <- runif(nrow(test$map)) str(test$map) # Reading again rds <- test$genotypes$rds test2 <- snp_attach(rds) str(test2$map) # new slot wasn't saved # Save it snp_save(test) # Reading again test3 <- snp_attach(rds) str(test3$map) # it is saved now ``` -------------------------------- ### Example: Compute and Inspect Correlation Matrix Source: https://privefl.github.io/bigsnpr/reference/snp_cor.html Demonstrates how to attach external data, compute a correlation matrix for the first 1000 SNPs, and display a subset of the results. ```R test <- snp_attachExtdata() G <- test$genotypes corr <- snp_cor(G, ind.col = 1:1000) corr[1:10, 1:10] ``` -------------------------------- ### Example Usage of bed_prodVec Source: https://privefl.github.io/bigsnpr/reference/bed_prodVec.html Demonstrates how to use the bed_prodVec function with a bed object and a vector. Missing values are handled by default. ```R bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") obj.bed <- bed(bedfile) y.col <- rep(1, ncol(obj.bed)) str(bed_prodVec(obj.bed, y.col)) #> num [1:517] 2841 2916 2842 2835 2872 ... ``` -------------------------------- ### SNP LD Split Example Source: https://privefl.github.io/bigsnpr/reference/snp_ldsplit.html Demonstrates how to use snp_ldsplit to split variants into LD blocks, visualize the trade-offs, and select an optimal split. ```R if (FALSE) { corr <- readRDS(url("https://www.dropbox.com/s/65u96jf7y32j2mj/spMat.rds?raw=1")) # adjust `THR_R2` depending on sample size used to compute corr # use e.g. 0.05 for small sample sizes, and 0.01 for large sample sizes THR_R2 <- 0.02 m <- ncol(corr) (SEQ <- round(seq_log(m / 30, m / 5, length.out = 10))) # replace `min_size` by e.g. 100 for larger data (res <- snp_ldsplit(corr, thr_r2 = THR_R2, min_size = 10, max_size = SEQ)) # add the variant block IDs corresponding to each split res$block_num <- lapply(res$all_size, function(.) rep(seq_along(.), .)) library(ggplot2) # trade-off cost / number of blocks qplot(n_block, cost, color = factor(max_size, SEQ), data = res) + theme_bw(14) + scale_y_log10() + theme(legend.position = "top") + labs(x = "Number of blocks", color = "Maximum block size", y = "Sum of squared correlations outside blocks") # trade-off cost / number of non-zero values qplot(perc_kept, cost, color = factor(max_size, SEQ), data = res) + theme_bw(14) + # scale_y_log10() + theme(legend.position = "top") + labs(x = "Percentage of non-zero values kept", color = "Maximum block size", y = "Sum of squared correlations outside blocks") # trade-off cost / sum of squared sizes qplot(cost2, cost, color = factor(max_size, SEQ), data = res) + theme_bw(14) + scale_y_log10() + geom_vline(xintercept = 0)+ theme(legend.position = "top") + labs(x = "Sum of squared blocks", color = "Maximum block size", y = "Sum of squared correlations outside blocks") ## Pick one solution and visualize blocks library(dplyr) all_ind <- res %>% arrange(cost2 * sqrt(5 + cost)) %>% print() %>% slice(1) %>% pull(all_last) ## 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, alpha = r2)) + theme_minimal() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), strip.background = element_blank(), strip.text.x = element_blank()) + scale_alpha_continuous(range = 0:1) + scale_x_continuous(expand = c(0.02, 0.02), minor_breaks = NULL, breaks = head(all_ind[[1]], -1) + 0.5) + facet_wrap(~ cut(i + y, 4), scales = "free", ncol = 1) + labs(x = "Position", y = NULL) } ``` -------------------------------- ### Example usage of same_ref Source: https://privefl.github.io/bigsnpr/reference/same_ref.html Demonstrates how to use the same_ref function with sample allele data. Missing values in the output can occur due to missing inputs or ambiguous allele matching. ```R same_ref(ref1 = c("A", "C", "T", "G", NA), alt1 = c("C", "T", "C", "A", "A"), ref2 = c("A", "C", "A", "A", "C"), alt2 = c("C", "G", "G", "G", "A")) #> [1] TRUE NA TRUE FALSE NA ``` -------------------------------- ### Example Usage of snp_plinkQC Source: https://privefl.github.io/bigsnpr/reference/snp_plinkQC.html Demonstrates how to use the snp_plinkQC function with custom QC parameters and output file naming. Ensure PLINK is downloaded and accessible via plink.path. ```R if (FALSE) { bedfile <- system.file("extdata", "example.bed", package = "bigsnpr") prefix <- sub_bed(bedfile) plink <- download_plink() test <- snp_plinkQC(plink.path = plink, prefix.in = prefix, prefix.out = tempfile(), file.type = "--bfile", # the default (for ".bed") maf = 0.05, geno = 0.05, mind = 0.05, hwe = 1e-10, autosome.only = TRUE) test } ``` -------------------------------- ### Download and Run PLINK Source: https://privefl.github.io/bigsnpr-extdoc/preprocessing.html Downloads the latest stable version of PLINK 1.9 and executes a command to check its version. Ensure PLINK is installed or use the download function. ```r plink <- download_plink("tmp-data") system(glue::glue( "{plink} --version" )) ``` ```text #> PLINK v1.9.0-b.7.7 64-bit (22 Oct 2024) ``` -------------------------------- ### Get FBM Backing File Path Source: https://privefl.github.io/bigsnpr-extdoc/working-with-an-fbm.html Retrieves the file path where the numeric data of the FBM is stored. ```R X$backingfile ``` -------------------------------- ### Example Usage of bed_scaleBinom Source: https://privefl.github.io/bigsnpr/reference/bed_scaleBinom.html Demonstrates how to use bed_scaleBinom to scale a bed object and how the output can be used with bed_randomSVD. The output of bed_scaleBinom is a data frame containing center and scale values. ```R bedfile <- system.file("extdata", "example-missing.bed", package = "bigsnpr") obj.bed <- bed(bedfile) str(bed_scaleBinom(obj.bed)) ``` ```R str(bed_randomSVD(obj.bed, bed_scaleBinom)) ``` -------------------------------- ### Example: LD clumping prioritizing higher MAF Source: https://privefl.github.io/bigsnpr/reference/snp_clumping.html Demonstrates how to use `snp_clumping` to select SNPs, prioritizing those with higher minor allele frequencies (MAFs). ```R test <- snp_attachExtdata() G <- test$genotypes # clumping (prioritizing higher MAF) ind.keep <- snp_clumping(G, infos.chr = test$map$chromosome, infos.pos = test$map$physical.pos, thr.r2 = 0.1) # keep most of them -> not much LD in this simulated dataset length(ind.keep) / ncol(G) ```