### Install SNPRelate Development Version from GitHub

Source: https://github.com/zhengxwen/snprelate/blob/master/README.md

Install the development version of SNPRelate and its dependency gdsfmt from GitHub. This method requires building from source and may need manual installation of dependencies.

```R
library("devtools")
install_github("zhengxwen/gdsfmt")
install_github("zhengxwen/SNPRelate")
```

--------------------------------

### Install SNPRelate from Bioconductor

Source: https://github.com/zhengxwen/snprelate/blob/master/README.md

Install the SNPRelate package from the Bioconductor repository using BiocManager. Ensure BiocManager is installed first.

```R
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("SNPRelate")
```

--------------------------------

### Load SNPRelate and gdsfmt Packages

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Load the required R packages for genomic data manipulation. Ensure these packages are installed before running.

```R
# Load the R packages: gdsfmt and SNPRelate
library(gdsfmt)
library(SNPRelate)
```

--------------------------------

### Summarize GDS File

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Use snpgdsSummary to get a summary of the created GDS file, including the number of samples and SNPs, and genotype storage mode.

```R
> # Summary
> snpgdsSummary("test.gds")
The file name: /Users/sts/Documents/Codes/Rpackages/SNPRelate_Sweave/output/test.gds 
The total number of samples: 279 
The total number of SNPs: 5000 
SNP genotypes are stored in individual-major mode (SNP X Sample).
```

--------------------------------

### Open GDS File and Get Population Information

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Opens a GDS file and reads population group information for samples. This is a prerequisite for many downstream analyses.

```R
> # Open the GDS file
> genofile <- snpgdsOpen(snpgdsExampleFileName())

> # Get population information
> #   or pop_code <- scan("pop.txt", what=character())
> #   if it is stored in a text file "pop.txt"
> pop_code <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"))
> # Display the first six values
> head(pop_code)
[1] "YRI" "YRI" "YRI" "YRI" "CEU" "CEU"
```

--------------------------------

### Read Genotype Data (First 5 Samples, First 3 SNPs)

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Reads a subset of genotype data from a GDS file, specifying the start and count for samples and SNPs. The output shows the genotype values (0, 1, 2, or NA).

```R
> g <- read.gdsn(index.gdsn(genofile, "genotype"), start=c(1,1), count=c(5,3))
     [,1] [,2] [,3]
[1,]    2    1    0
[2,]    1    1    0
[3,]    2    1    1
[4,]    2    1    1
[5,]    0    0    0
```

--------------------------------

### Get Genotype Data by Sample and SNP IDs

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Retrieves genotype data using specific sample and SNP IDs. This function replaces a value of 3 with NA.

```R
> g <- snpgdsGetGeno(genofile, sample.id=..., snp.id=...)
```

--------------------------------

### Get Genotype Attribute

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Retrieves attributes of the genotype data, such as sample or SNP order, which indicates the data's dimension mode (individual-major or SNP-major).

```R
> get.attr.gdsn(index.gdsn(genofile, "genotype"))
$sample.order
NULL
```

--------------------------------

### Prepare Sample Loadings Data with Population Info

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Creates a data frame that includes sample IDs, population codes, and the first two eigenvector values. This is necessary for coloring plots by population group.

```R
> # Get sample id
> sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
> # Get population information
> #   or pop_code <- scan("pop.txt", what=character())
> #   if it is stored in a text file "pop.txt"
> pop_code <- read.gdsn(index.gdsn(genofile, "sample.annot/pop.group"))
> # assume the order of sample IDs is as the same as population codes
> head(cbind(sample.id, pop_code))
     sample.id pop_code
[1,] "NA19152" "YRI"   
[2,] "NA19139" "YRI"   
[3,] "NA18912" "YRI"   
[4,] "NA19160" "YRI"   
[5,] "NA07034" "CEU"   
[6,] "NA07055" "CEU"   

> # Make a data.frame
> tab <- data.frame(sample.id = pca$sample.id,
    pop = factor(pop_code)[match(pca$sample.id, sample.id)],
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
> head(tab)
  sample.id pop         EV1         EV2
1   NA19152 YRI  0.08411287  0.01226860
2   NA19139 YRI  0.08360644  0.01085849
3   NA18912 YRI  0.08110808  0.01184524
4   NA19160 YRI  0.08680864  0.01447106
5   NA07034 CEU -0.03109761 -0.07709255
6   NA07055 CEU -0.03228450 -0.08155730
```

--------------------------------

### Advanced GDS File Creation with gdsfmt

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

This snippet demonstrates advanced GDS file creation and manipulation using functions from the gdsfmt package. It covers creating a new file, adding attributes, adding variables (like sample IDs, SNP IDs, positions, chromosomes, and alleles), and writing genotype data in either SNP-by-sample or sample-by-SNP format. It also shows how to add sample and SNP annotations and set chromosome code attributes.

```R
# Create a new GDS file
newfile <- createfn.gds("your_gds_file.gds")

# add a flag
put.attr.gdsn(newfile$root, "FileFormat", "SNP_ARRAY")

# Add variables
add.gdsn(newfile, "sample.id", sample.id)
add.gdsn(newfile, "snp.id", snp.id)
add.gdsn(newfile, "snp.position", snp.position)
add.gdsn(newfile, "snp.chromosome", snp.chromosome)
add.gdsn(newfile, "snp.allele", c("A/G", "T/C", ...))

#####################################################################
# Create a snp-by-sample genotype matrix

# Add genotypes
var.geno <- add.gdsn(newfile, "genotype",
    valdim=c(length(snp.id), length(sample.id)), storage="bit2")

# Indicate the SNP matrix is snp-by-sample
put.attr.gdsn(var.geno, "snp.order")

# Write SNPs into the file sample by sample
for (i in 1:length(sample.id))
{
    g <- ...
    write.gdsn(var.geno, g, start=c(1,i), count=c(-1,1))
}


#####################################################################
# OR, create a sample-by-snp genotype matrix

# Add genotypes
var.geno <- add.gdsn(newfile, "genotype",
    valdim=c(length(sample.id), length(snp.id)), storage="bit2")

# Indicate the SNP matrix is sample-by-snp
put.attr.gdsn(var.geno, "sample.order")

# Write SNPs into the file sample by sample
for (i in 1:length(snp.id))
{
    g <- ...
    write.gdsn(var.geno, g, start=c(1,i), count=c(-1,1))
}


# Get a description of chromosome codes
#   allowing to define a new chromosome code, e.g., snpgdsOption(Z=27)
option <- snpgdsOption()
var.chr <- index.gdsn(newfile, "snp.chromosome")
put.attr.gdsn(var.chr, "autosome.start", option$autosome.start)
put.attr.gdsn(var.chr, "autosome.end", option$autosome.end)
for (i in 1:length(option$chromosome.code))
{
    put.attr.gdsn(var.chr, names(option$chromosome.code)[i],
        option$chromosome.code[[i]])
}

# Add your sample annotation
samp.annot <- data.frame(sex = c("male", "male", "female", ...),
    pop.group = c("CEU", "CEU", "JPT", ...), ...)
add.gdsn(newfile, "sample.annot", samp.annot)

# Add your SNP annotation
snp.annot <- data.frame(pass=c(TRUE, TRUE, FALSE, FALSE, TRUE, ...), ...)
add.gdsn(newfile, "snp.annot", snp.annot)


# Close the GDS file
closefn.gds(newfile)
```

--------------------------------

### Configure Makevars for Intel SIMD Intrinsics

Source: https://github.com/zhengxwen/snprelate/blob/master/README.md

Set CFLAGS and CXXFLAGS in ~/.R/Makevars to enable native CPU optimizations for C and C++ code when compiling R packages.

```sh
## for C code
CFLAGS=-g -O3 -march=native -mtune=native
## for C++ code
CXXFLAGS=-g -O3 -march=native -mtune=native
```

--------------------------------

### Create Data Frame for Sample Loadings (No Population Info)

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Creates a data frame containing sample IDs and their first two eigenvector values (principal components). This is useful for plotting sample positions in a reduced dimensional space.

```R
> # make a data.frame
> tab <- data.frame(sample.id = pca$sample.id,
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
> head(tab)
  sample.id         EV1         EV2
1   NA19152  0.08411287  0.01226860
2   NA19139  0.08360644  0.01085849
3   NA18912  0.08110808  0.01184524
4   NA19160  0.08680864  0.01447106
5   NA07034 -0.03109761 -0.07709255
6   NA07055 -0.03228450 -0.08155730
```

--------------------------------

### Run Principal Component Analysis

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Calculates the principal components for genotype data. Ensure the genofile object is properly loaded before running.

```R
> # Run PCA
> pca <- snpgdsPCA(genofile)
Principal Component Analysis (PCA) on SNP genotypes:
Removing 365 non-autosomal SNPs.
Removing 1 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 279 samples, 8722 SNPs
	Using 1 (CPU) core
PCA:	the sum of all working genotypes (0, 1 and 2) = 2446510
PCA:	Thu Jul 17 17:48:04 2014	0%
PCA:	Thu Jul 17 17:48:04 2014	100%
PCA:	Thu Jul 17 17:48:04 2014	Begin (eigenvalues and eigenvectors)
PCA:	Thu Jul 17 17:48:04 2014	End (eigenvalues and eigenvectors)
```

--------------------------------

### Plot Sample Loadings (No Population Info)

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the first eigenvector against the second eigenvector for all samples. This visualization helps identify population structure or outliers.

```R
> # Draw
> plot(tab$EV2, tab$EV1, xlab="eigenvector 2", ylab="eigenvector 1")
```

--------------------------------

### Create GDS File with snpgdsCreateGeno

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Use snpgdsCreateGeno to create a GDS file from a numeric genotype matrix. Ensure sample IDs, SNP IDs, chromosome, position, and allele information are provided. The snpfirstdim argument controls the matrix orientation.

```R
> # Load data
> data(hapmap_geno)
> # Create a gds file
> snpgdsCreateGeno("test.gds", genmat = hapmap_geno$genotype,
    sample.id = hapmap_geno$sample.id, snp.id = hapmap_geno$snp.id,
    snp.chromosome = hapmap_geno$snp.chromosome,
    snp.position = hapmap_geno$snp.position,
    snp.allele = hapmap_geno$snp.allele, snpfirstdim=TRUE)

> # Open the GDS file
> (genofile <- snpgdsOpen("test.gds"))
File: /Users/sts/Documents/Codes/Rpackages/SNPRelate_Sweave/output/test.gds
+ 	[  ] *
|--+ sample.id	{ VStr8 279 ZIP(29.89%) }
|--+ snp.id	{ VStr8 1000 ZIP(42.42%) }
|--+ snp.position	{ Float64 1000 ZIP(55.97%) }
|--+ snp.chromosome	{ Int32 1000 ZIP(2.00%) }
|--+ snp.allele	{ VStr8 1000 ZIP(13.85%) }
|--+ genotype	{ Bit2 1000x279 } *

> # Close the GDS file
> snpgdsClose(genofile)
```

--------------------------------

### Summarize VCF-converted GDS File

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Use snpgdsSummary to review the GDS file created from a VCF file. This shows the number of samples, SNPs, and the genotype storage mode.

```R
> # Summary
> snpgdsSummary("test.gds")
The file name: /Users/sts/Documents/Codes/Rpackages/SNPRelate_Sweave/output/test.gds 
The total number of samples: 3 
The total number of SNPs: 2 
SNP genotypes are stored in SNP-major mode (Sample X SNP).
```

--------------------------------

### Open and Inspect a GDS File

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Open a GDS file for reading and display its structure, including sample IDs, SNP information, and genotype data. This helps understand the data organization within the file.

```R
> # Open a GDS file
> (genofile <- snpgdsOpen(snpgdsExampleFileName()))
File: /library/SNPRelate/extdata/hapmap_geno.gds
+ 	[  ] *
|--+ sample.id	{ FStr8 279 ZIP(23.10%) }
|--+ snp.id	{ Int32 9088 ZIP(34.76%) }
|--+ snp.rs.id	{ FStr8 9088 ZIP(42.66%) }
|--+ snp.position	{ Int32 9088 ZIP(94.73%) }
|--+ snp.chromosome	{ UInt8 9088 ZIP(0.94%) } *
|--+ snp.allele	{ FStr8 9088 ZIP(14.45%) }
|--+ genotype	{ Bit2 279x9088 } *
|--+ sample.annot	[ data.frame ] *
   |--+ sample.id	{ FStr8 279 ZIP(23.10%) }
   |--+ family.id	{ FStr8 279 ZIP(28.37%) }
   |--+ father.id	{ FStr8 279 ZIP(12.98%) }
   |--+ mother.id	{ FStr8 279 ZIP(12.86%) }
   |--+ sex	{ FStr8 279 ZIP(28.32%) }
   |--+ pop.group	{ FStr8 279 ZIP(7.89%) }
```

--------------------------------

### Plot Sample Loadings with Population Colors

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the first eigenvector against the second, coloring points by population group. A legend is included to identify the population codes.

```R
> # Draw
> plot(tab$EV2, tab$EV1, col=as.integer(tab$pop),
    xlab="eigenvector 2", ylab="eigenvector 1")
> legend("topleft", legend=levels(tab$pop), pch="o", col=1:nlevels(tab$pop))
```

--------------------------------

### Convert PLINK Binary to GDS

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Use snpgdsBED2GDS to convert PLINK binary files to GDS format. Specify the paths to your .bed, .fam, and .bim files, and the output GDS file name.

```R
> # The PLINK BED file, using the example in the SNPRelate package
> bed.fn <- system.file("extdata", "plinkhapmap.bed", package="SNPRelate")
> fam.fn <- system.file("extdata", "plinkhapmap.fam", package="SNPRelate")
> bim.fn <- system.file("extdata", "plinkhapmap.bim", package="SNPRelate")
```

```R
> bed.fn <- "C:/your_folder/your_plink_file.bed"
> fam.fn <- "C:/your_folder/your_plink_file.fam"
> bim.fn <- "C:/your_folder/your_plink_file.bim"
> # Convert
> snpgdsBED2GDS(bed.fn, fam.fn, bim.fn, "test.gds")
Start snpgdsBED2GDS ...
	BED file: "/library/SNPRelate/extdata/plinkhapmap.bed" in the individual-major mode (SNP X Sample)
	FAM file: "/library/SNPRelate/extdata/plinkhapmap.fam", DONE.
	BIM file: "/library/SNPRelate/extdata/plinkhapmap.bim", DONE.
Thu Jul 17 17:48:04 2014 	store sample id, snp id, position, and chromosome.
	start writing: 279 samples, 5000 SNPs ...
 	Thu Jul 17 17:48:04 2014	0%
 	Thu Jul 17 17:48:04 2014	100%
Thu Jul 17 17:48:04 2014 	Done.
```

--------------------------------

### Visualize IBS Matrix with Heat Map

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Visualize the calculated IBS matrix using a heat map with the lattice package. The order of individuals can be adjusted based on population codes for better visualization.

```R
> library(lattice)
> L <- order(pop_code)
> levelplot(ibs$ibs[L, L], col.regions = terrain.colors)
```

--------------------------------

### Calculate Percent Variation Explained by PCs

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Calculates and displays the percentage of total variation accounted for by the first 16 principal components. This helps in understanding the explanatory power of each component.

```R
> pc.percent <- 100 * pca$eigenval[1:16]/sum(pca$eigenval)
> pc.percent
 [1] 12.2251756  5.8397767  1.0111423  0.9475002  0.8437720  0.7363588
 [7]  0.7350545  0.7243414  0.7200081  0.7047742  0.7033253  0.6943755
[13]  0.6847090  0.6780401  0.6720635  0.6692861
```

--------------------------------

### Plot IBD Coefficients (MoM)

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the k0 (proportion of homozygous by descent) against k1 (proportion of heterozygous by descent) values obtained from the MoM IBD estimation.

```R
> plot(ibd.coeff$k0, ibd.coeff$k1, xlim=c(0,1), ylim=c(0,1),
    xlab="k0", ylab="k1", main="YRI samples (MoM)")
> lines(c(0,1), c(1,0), col="red", lty=2)
```

--------------------------------

### Summarize GDS File Information

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Display summary statistics of a GDS file, including sample and SNP counts, and genotype storage mode. This is useful for initial data inspection.

```R
> snpgdsSummary(snpgdsExampleFileName())
The file name: /library/SNPRelate/extdata/hapmap_geno.gds 
The total number of samples: 279 
The total number of SNPs: 9088 
SNP genotypes are stored in SNP-major mode (Sample X SNP).
```

--------------------------------

### Determine Groups by Population Information using Hierarchical Clustering

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Determine groups of individuals using hierarchical clustering based on provided population information. The dendrogram is plotted, and a legend is added to indicate population groups.

```R
> # Determine groups of individuals by population information
> rv2 <- snpgdsCutTree(ibs.hc, samp.group=as.factor(pop_code))
Create 4 groups.

> plot(rv2$dendrogram, leaflab="none", main="HapMap Phase II")
> legend("topright", legend=levels(race), col=1:nlevels(race), pch=19, ncol=4)
```

--------------------------------

### Read Sample Population Information

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Reads population group information for samples from the GDS file and displays a frequency table of the populations.

```R
> # Read population information
> pop <- read.gdsn(index.gdsn(genofile, path="sample.annot/pop.group"))
> table(pop)
pop
CEU HCB JPT YRI 
 92  47  47  93 
```

--------------------------------

### Plot Pairs of Principal Components

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Generates a scatter plot matrix for the first four principal components, with points colored by population. Labels indicate the PC number and the percentage of variance explained.

```R
> lbls <- paste("PC", 1:4, "\n", format(pc.percent[1:4], digits=2), "%", sep="")
> pairs(pca$eigenvect[,1:4], col=tab$pop, labels=lbls)
```

--------------------------------

### Read Sample IDs from GDS File

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Reads sample IDs from a GDS file and filters them based on population code.

```R
> # YRI samples
> sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
> YRI.id <- sample.id[pop_code == "YRI"]
```

--------------------------------

### Plot IBD Coefficients (MLE)

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the k0 against k1 values obtained from the MLE IBD estimation, similar to the MoM plot.

```R
> plot(ibd.coeff$k0, ibd.coeff$k1, xlim=c(0,1), ylim=c(0,1),
    xlab="k0", ylab="k1", main="YRI samples (MLE)")
> lines(c(0,1), c(1,0), col="red", lty=2)
```

--------------------------------

### Plot IBS0 vs Kinship Coefficient (KING-robust)

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the Proportion of Zero IBS (IBS0) against the Estimated Kinship Coefficient, derived from the KING-robust method. This visualization helps in understanding pairwise relationships.

```R
> plot(dat$IBS0, dat$kinship, xlab="Proportion of Zero IBS",
    ylab="Estimated Kinship Coefficient (KING-robust)")
```

--------------------------------

### Infer Relationships using KING Method

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Infers within- and between-family relationships using the KING-robust method, which is effective in the presence of population stratification. It requires sample IDs and family IDs.

```R
> # Incorporate with pedigree information
> family.id <- read.gdsn(index.gdsn(genofile, "sample.annot/family.id"))
> family.id <- family.id[match(YRI.id, sample.id)]
> table(family.id)
family.id
101 105 112 117  12  13  16  17  18  23  24  28   4  40  42  43  45  47  48   5 
  3   3   3   4   4   3   3   3   3   3   4   3   3   3   3   3   3   3   3   3 
 50  51  56  58  60  71  72  74  77   9 
  3   3   3   3   3   3   3   3   3   3 
```

```R
> ibd.robust <- snpgdsIBDKING(genofile, sample.id=YRI.id, family.id=family.id)
IBD analysis (KING method of moment) on SNP genotypes:
Removing 365 non-autosomal SNPs.
Removing 563 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 93 samples, 8160 SNPs
	Using 1 (CPU) core
# of families: 30, and within- and between-family relationship are estimated differently.
Working space: 93 samples, 8160 SNPs
	Using 1 CPU core.
Relationship inference in the presence of population stratification.
KING IBD:	the sum of all working genotypes (0, 1 and 2) = 755648
KING IBD:	Thu Jul 17 17:48:26 2014	0%
KING IBD:	Thu Jul 17 17:48:26 2014	100%
```

```R
> names(ibd.robust)
[1] "sample.id" "snp.id"    "afreq"     "IBS0"      "kinship"  
```

```R
> # Pairs of individuals
> dat <- snpgdsIBDSelection(ibd.robust)
> head(dat)
      ID1     ID2       IBS0      kinship
1 NA19152 NA19139 0.05504926 -0.005516960
2 NA19152 NA18912 0.05738916 -0.003658537
3 NA19152 NA19160 0.06230760 -0.034086156
4 NA19152 NA18515 0.05602758  0.007874016
5 NA19152 NA19222 0.05923645 -0.012668574
6 NA19152 NA18508 0.05561722  0.002216848
```

--------------------------------

### Estimate IBD Coefficients using MLE

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Estimates IBD coefficients using Maximum Likelihood Estimation (MLE). A random subset of SNPs can be selected for this analysis. The output can be converted to a data frame.

```R
> # Estimate IBD coefficients
> set.seed(1000)
> snp.id <- sample(snpset.id, 5000)  # random 5000 SNPs
> ibd <- snpgdsIBDMLE(genofile, sample.id=YRI.id, snp.id=snp.id,
    maf=0.05, missing.rate=0.05)
```

```R
> # Make a data.frame
> ibd.coeff <- snpgdsIBDSelection(ibd)
```

--------------------------------

### Perform Hierarchical Clustering on IBS Distances and Determine Groups

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Perform hierarchical clustering on the IBS distance matrix and automatically determine groups of individuals using a permutation score. The resulting dendrogram and group assignments can be visualized and tabulated.

```R
> set.seed(100)
> ibs.hc <- snpgdsHCluster(snpgdsIBS(genofile, num.thread=2))
Identity-By-State (IBS) analysis on SNP genotypes:
Removing 365 non-autosomal SNPs.
Removing 1 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 279 samples, 8722 SNPs
	Using 2 (CPU) cores
IBS:	the sum of all working genotypes (0, 1 and 2) = 2446510
IBS:	Thu Jul 17 17:48:28 2014	0%
IBS:	Thu Jul 17 17:48:28 2014	100%

> # Determine groups of individuals automatically
> rv <- snpgdsCutTree(ibs.hc)
Determine groups by permutation (Z threshold: 15, outlier threshold: 5):
Create 3 groups.

> plot(rv$dendrogram, leaflab="none", main="HapMap Phase II")
> table(rv$samp.group)
```

--------------------------------

### Read SNP IDs

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Extracts and displays the first few SNP IDs from the GDS file.

```R
> # Take out snp.id
> head(read.gdsn(index.gdsn(genofile, "snp.id")))
[1] 1 2 3 4 5 6
```

--------------------------------

### Calculate Genome-Wide Average IBS Pairwise Identities

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Use snpgdsIBS to compute an n-by-n matrix of genome-wide average IBS pairwise identities for a given set of SNP genotypes. This function supports multithreading for faster computation.

```R
> ibs <- snpgdsIBS(genofile, num.thread=2)
Identity-By-State (IBS) analysis on SNP genotypes:
Removing 365 non-autosomal SNPs.
Removing 1 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 279 samples, 8722 SNPs
	Using 2 (CPU) cores
IBS:	the sum of all working genotypes (0, 1 and 2) = 2446510
IBS:	Thu Jul 17 17:48:26 2014	0%
IBS:	Thu Jul 17 17:48:26 2014	100%
```

--------------------------------

### Estimate IBD Coefficients using PLINK MoM

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Estimates IBD coefficients using the PLINK method of moments. It filters SNPs based on MAF and missing rate before computation. The output can be converted to a data frame for further analysis.

```R
> # Estimate IBD coefficients
> ibd <- snpgdsIBDMoM(genofile, sample.id=YRI.id, snp.id=snpset.id,
    maf=0.05, missing.rate=0.05)
IBD analysis (PLINK method of moment) on SNP genotypes:
Removing 1285 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 93 samples, 5262 SNPs
	Using 1 (CPU) core
PLINK IBD:	the sum of all working genotypes (0, 1 and 2) = 484520
PLINK IBD:	Thu Jul 17 17:48:26 2014	0%
PLINK IBD:	Thu Jul 17 17:48:26 2014	100%
```

```R
> # Make a data.frame
> ibd.coeff <- snpgdsIBDSelection(ibd)
> head(ibd.coeff)
      ID1     ID2        k0         k1     kinship
1 NA19152 NA19139 0.9548539 0.04514610 0.011286524
2 NA19152 NA18912 1.0000000 0.00000000 0.000000000
3 NA19152 NA19160 1.0000000 0.00000000 0.000000000
4 NA19152 NA18515 0.9234541 0.07654590 0.019136475
5 NA19152 NA19222 1.0000000 0.00000000 0.000000000
6 NA19152 NA18508 0.9833803 0.01661969 0.004154922
```

--------------------------------

### Read SNP RS IDs

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Extracts and displays the first few reference SNP IDs (rs IDs) from the GDS file.

```R
> # Take out snp.rs.id
> head(read.gdsn(index.gdsn(genofile, "snp.rs.id")))
[1] "rs1695824"  "rs13328662" "rs4654497"  "rs10915489" "rs12132314"
[6] "rs12042555"
```

--------------------------------

### Perform Multidimensional Scaling Analysis on IBS Distances

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Conduct multidimensional scaling (MDS) analysis on the IBS distance matrix to visualize relationships between individuals in a 2D space. The plot can be colored by population groups and includes a legend.

```R
> loc <- cmdscale(1 - ibs$ibs, k = 2)
> x <- loc[, 1]; y <- loc[, 2]
> race <- as.factor(pop_code)

> plot(x, y, col=race, xlab = "", ylab = "",
    main = "Multidimensional Scaling Analysis (IBS Distance)")
> legend("topleft", legend=levels(race), text.col=1:nlevels(race))
```

--------------------------------

### Convert VCF to GDS

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Use snpgdsVCF2GDS to convert VCF files to GDS format. The 'method' argument can be set to 'biallelic.only' to process only biallelic SNPs.

```R
> # The VCF file, using the example in the SNPRelate package
> vcf.fn <- system.file("extdata", "sequence.vcf", package="SNPRelate")
```

```R
> vcf.fn <- "C:/your_folder/your_vcf_file.vcf"
> # Reformat
> snpgdsVCF2GDS(vcf.fn, "test.gds", method="biallelic.only")
VCF format ---> GDS SNP format:
Parsing "/library/SNPRelate/extdata/sequence.vcf" ...
	Import 2 variants.
Optimize the access efficiency ...
Clean up the fragments of GDS file:
	open the file "test.gds" (size: 2011).
	# of fragments in total: 22.
	save it to "test.gds.tmp".
	rename "test.gds.tmp" (size: 1987).
	# of fragments in total: 20.
```

--------------------------------

### Close GDS File

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Closes the connection to the GDS file after data has been read.

```R
> # Close the GDS file
> snpgdsClose(genofile)
```

--------------------------------

### Retrieve Chromosome Coding Attributes

Source: https://github.com/zhengxwen/snprelate/wiki/Preparing-Data

Access and display the attributes associated with the 'snp.chromosome' variable in a GDS file, which define the numeric coding for chromosomes. This is useful for understanding how chromosomes are represented.

```R
> # Get the attributes of chromosome coding
> get.attr.gdsn(index.gdsn(genofile, "snp.chromosome"))
$autosome.start
[1] 1
$autosome.end
[1] 22
$X
[1] 23
$XY
[1] 24
$Y
[1] 25
$M
[1] 26
$MT
[1] 26
```

--------------------------------

### LD-based SNP Pruning

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Performs SNP pruning based on linkage disequilibrium (LD) to select a subset of SNPs in approximate LD. This is useful for principal component analysis and relatedness analysis.

```R
> set.seed(1000)
> # Try different LD thresholds for sensitivity analysis
> snpset <- snpgdsLDpruning(genofile, ld.threshold=0.2)
SNP pruning based on LD:
Removing 365 non-autosomal SNPs.
Removing 1 SNPs (monomorphic, < MAF, or > missing rate)
Working space: 279 samples, 8722 SNPs
	Using 1 (CPU) core
	Sliding window: 500000 basepairs, Inf SNPs
	|LD| threshold: 0.2
Chromosome 1: 75.42%, 540/716
Chromosome 2: 72.24%, 536/742
Chromosome 3: 74.71%, 455/609
Chromosome 4: 73.31%, 412/562
Chromosome 5: 77.03%, 436/566
Chromosome 6: 75.58%, 427/565
Chromosome 7: 75.42%, 356/472
Chromosome 8: 71.31%, 348/488
Chromosome 9: 77.88%, 324/416
Chromosome 10: 74.33%, 359/483
Chromosome 11: 77.40%, 346/447
Chromosome 12: 76.81%, 328/427
Chromosome 13: 75.58%, 260/344
Chromosome 14: 76.95%, 217/282
Chromosome 15: 76.34%, 200/262
Chromosome 16: 72.66%, 202/278
Chromosome 17: 74.40%, 154/207
Chromosome 18: 73.68%, 196/266
Chromosome 19: 85.00%, 102/120
Chromosome 20: 71.62%, 164/229
Chromosome 21: 76.98%, 97/126
Chromosome 22: 75.86%, 88/116
6547 SNPs are selected in total.

> names(snpset)
 [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9" 
[10] "chr10" "chr11" "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18"
[19] "chr19" "chr20" "chr21" "chr22" 
> head(snpset$chr1)  # snp.id
[1]  1  2  4  5  7 10
> # Get all selected snp id
> snpset.id <- unlist(snpset)

```

--------------------------------

### Calculate SNP Correlations with Eigenvectors

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Calculates the correlation between SNP genotypes and the top principal components. This helps identify which SNPs contribute most to the observed variation.

```R
> # Get chromosome index
> chr <- read.gdsn(index.gdsn(genofile, "snp.chromosome"))
> CORR <- snpgdsPCACorr(pca, genofile, eig.which=1:4)
SNP correlations:
Working space: 279 samples, 9088 SNPs
	Using 1 CPU core.
	Using the top 32 eigenvectors.
SNP Correlation:	the sum of all working genotypes (0, 1 and 2) = 2553065
SNP Correlation:	Thu Jul 17 17:48:05 2014	0%
SNP Correlation:	Thu Jul 17 17:48:05 2014	100%
```

--------------------------------

### Plot SNP Correlations

Source: https://github.com/zhengxwen/snprelate/wiki/Data-Analysis

Plots the absolute correlation values between SNPs and the first three principal components. SNPs are colored by chromosome, and the y-axis represents the correlation strength.

```R
> par( mfrow=c(3,1))
> for (i in 1:3)
  {
      plot(abs(CORR$snpcorr[i,]), ylim=c(0,1), xlab="SNP Index",
          ylab=paste("PC", i), col=chr, pch="+`)
  }
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.