2.4.1.1 Examples

clq1 = CLQD(geno[, 1:100], SNPinfo[1:100,])

2.4.1.2 Values

CLQD returns a vector of bin numbers of SNPs. NA means that the corresponding SNP is a singleton, i.e. that SNP is an LD bin of size 1.

# CLQ results for 300 SNPs
head(clq1, n = 20)

##  [1] 10 NA NA 14 12 NA 15 NA NA  8 NA 10 NA 10 NA 11 NA NA NA  4

table(clq1)

## clq1
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 
##  3 10  2  2  3  3  5  2  7 11  2  2  2  2  2

# n of singletons
sum(is.na(clq1))

## [1] 42

2.4.1.3 Details

CLQD(geno[, 1:100], SNPinfo[1:100,], LD = "Dprime")

The algorithm first constructs a graph in which the vertices are SNPs and each edge is an unordered pair of two SNPs. Each edge weight is given by LD value between two SNPs. CLQD uses \(r^2\) or D’ as the LD measurement. If you want to measure \(r^2\), set LD = "r2"(default), otherwise set LD = "Dprime". When LD = "r2", we have to specify a threshold for \(|r|\) by using the parameter CLQcut. The default is CLQcut = 0.5. For the D’, we define valid edges following definition of ‘strong LD’ suggested by Gabriel et al., (2002).

CLQD(geno[, 1:100], SNPinfo[1:100,], CLQcut = 0.7, hrstType = "fast", 
     CLQmode = "maximal")

CLQ algorithm first detects maximal cliques in the constructed graph and then takes cliques which include SNPs not yet selected as LD bins greedily. There are two clique selection methods. If you want to give a priority to a clique of larger size, set as CLQmode = "maximal", or if you want to give a priority to a cluster of denser LD bins (that is, the average distance of SNPs in a bin is smaller), set as CLQmode = "density"(default).

For the worst case, memory and time to listing all maximal cliques could be exponentially increased. Although, we can find all maximal cliques in reasonable time and memory usage for the most graphs corresponding to the LD structures. However, there are SNP sequence regions in which the time and memory usage for detecting maximal cliques extremely increasing, and we suggested heuristic algorithm for these regions. There are two kind of heuristic algorithm, “fast” and “near-nonhrst”(default). You can choose which algorithm to be used by setting the option hrstType such as hrstType = "fast" or hrstType = "near-nonhrst". “fast” algorithm needs less memory usage and time than the “near-nonhrst” algorithm. The “near-nonhrst” algorithm needs more memory and time than “fast” algorithm, but it returns the results similar to results obtained by non-heuristic algorithm. Adopting the “near-nonhrst” algorithm should specify a heuristic parameter hrstParam that specifies how similar to the non-heuristic result. Entering a high value produces results similar to non-heuristic results. It is recommended to set at least 150, and the default value is 200. To use the non-heuristic algorithm, set as hrstType = "nonhrst".

If the distance of two consecutive SNPs in a clique is too far, we can divide the clique into two cliques. The parameter for the threshold is clstgap and the default is 40000(bp).

2.4.2.1 Examples

# Use the options `geno` and `SNPinfo` to input additive genotype data and SNP information data respectively.
# use r2 measure (default)
res1 = BigLD(geno = geno[,1:1000], SNPinfo = SNPinfo[1:1000,]) 

## working chromosome name: 21

## remove SNPs with MAF < 0.05

## Number of SNPs after prunning SNPs with MAF< 0.05 : 501

## start to split the whole sequence into sub-task regions

## there is only one sub-region!

## 2019-05-03 01:31:13 | chr21:16040232-16196907 | sub-region: 1/1

## CLQ done!
constructLDblock done!

## 
## BigLD done!

#use D' measure
res1_dp = BigLD(geno = geno[,1:1000], SNPinfo = SNPinfo[1:1000,]
                ,LD = "Dprime") 

2.4.2.2 Values

BigLD returns a data frame which contains 7 information of each block (chromosome name, indices of the first SNP and last SNP, rsID of the first SNP and last SNP, bp of the first SNP and last SNP)

head(res1)

##   chr start.index end.index start.rsID   end.rsID start.bp   end.bp
## 1  21           1        14 rs59504523  rs2822825 16040232 16042876
## 2  21          20        23  rs9808806 rs13052089 16043928 16044788
## 3  21          26        93  rs2822826  rs2822833 16045615 16057487
## 4  21          95       115   rs988831   rs960479 16057936 16060592
## 5  21         125       132   rs467116 rs12626958 16062725 16063744
## 6  21         133       147   rs458700   rs458705 16064006 16065336

2.4.2.3 Details

# When you have to input files directly, use the parameter `genofile` (and `SNPinfofile`) instead of `geno` and `SNPinfo`.
res2 = BigLD(genofile = "geno.vcf")
res3 = BigLD(genofile = "geno.ped", SNPinfofile = "geno.map")
# change LD measure, hrstParam
res4 = BigLD(geno = geno, SNPinfo = SNPinfo, LD = "Dprime", hrstParam = 150)

The input consists of genotype and SNP information. The function also can directly load the data from the files of PLINK format or VCF (.ped, .map, .traw, .raw, .vcf). If the file sizes (.ped, .traw, .raw, .vcf) are larger than 10Gb, we strongly recommand to use geno and SNPinfo options after reformatting your data instead of loading large files file directly. As with CLQD function, you can choose an LD measure to use (LD), r2 threshold (CLQcut), heuristic algorithm (hrstType, hrstParam), CLQ methods (CLQmode), threshold for distance of consecutive SNPs in a clique (clstgap).

res5 = BigLD(geno = geno, SNPinfo = SNPinfo, MAFcut = 0.1, appendRare = TRUE)

BigLD algorithm construct LD blocks using the common variants after filtering rare variants by the option MAFcut(default 0.05). To include a rare variant in the LD block results, add the appendRare = TRUE option. Rare variants with maf \(\leq\) MAFcut could be appended to a proper LD block already constructed by common variants.

We can force specified SNP to be the last SNP in a block using cutByForce option.

cutlist = rbind(c(21, "rs440600", 16051956), c(21, "rs9979041", 16055738))
res6 = BigLD(geno = geno[, 1:100], SNPinfo = SNPinfo[1:100,], cutByForce = cutlist)

print(cutlist)

##      [,1] [,2]        [,3]      
## [1,] "21" "rs440600"  "16051956"
## [2,] "21" "rs9979041" "16055738"

head(res6)

##   chr start.index end.index start.rsID   end.rsID start.bp   end.bp
## 1  21           1        93 rs59504523  rs2822833 16040232 16057487
## 2  21          95       100   rs988831 rs59861729 16057936 16058409

If you want to obtain BigLD results of some specific region of an input data, you can speicify the region to calculate by using the options chrN, startbp, and endbp.

res7 = BigLD(geno = geno, SNPinfo = SNPinfo, startbp = 16058400, endbp = 16076500)

res7

##   chr start.index end.index start.rsID   end.rsID start.bp   end.bp
## 1  21           1        16 rs59861729   rs960479 16058409 16060592
## 2  21          26        99   rs467116 rs13051994 16062725 16076420

2.4.3.1 Examples

# gene based GPART using the pre-calculated BigLD result 'res1'
# The input data of GPART must be the same as the input data used to obtain 'res1'
# note that the `res1` is obtained by using the first 1000 SNPs in geno.
# default minsize = 4, maxsize = 50
Gres0 = GPART(geno = geno[,1:1000], SNPinfo = SNPinfo[1:1000,], 
              geneinfo = geneinfo,BigLDresult = res1, 
              minsize = 4, maxsize = 50)

## Number of SNPs after prunning SNPs with MAF< 0.05 : 501

## Use the inputted BigLD result

## GPART algorithm: geneBased

## working chromosome name: 21

## chr 21 done!

2.4.3.2 Values

GPART() returns data frame which contains 9 information of each partition (chromosome, indices of the first SNP and last SNP in the inputted data, rsIDs of the first SNP and last SNP, basepair positions of the first SNP and last SNP, blocksize, Name of a block)

# results of gene-based GPART
head(Gres0)

##   chr start.index end.index start.rsID   end.rsID start.bp   end.bp
## 1  21           1         5 rs59504523 rs13052089 16040232 16044788
## 2  21           6        43  rs2822826  rs2822833 16045615 16057487
## 3  21          44        53   rs988831   rs960479 16057936 16060592
## 4  21          54        57   rs467116 rs12626958 16062725 16063744
## 5  21          58        68   rs458700   rs458705 16064006 16065336
## 6  21          69       116   rs460403  rs2822849 16066007 16078541
##   blocksize                              Name
## 1         5 inter-AF165138.7:AF127936.3-part1
## 2        38 inter-AF165138.7:AF127936.3-part2
## 3        10 inter-AF165138.7:AF127936.3-part3
## 4         4 inter-AF165138.7:AF127936.3-part4
## 5        11 inter-AF165138.7:AF127936.3-part5
## 6        48 inter-AF165138.7:AF127936.3-part6

< Naming Examples >

• before-Gene1-part1 : 1st block located before 1st gene (No overlapping LD block with the Gene1).
• before-Gene1-part2 : 2nd block located at before 1st gene region.
• Gene2/Gene3-part1 : 1st block overlapping Gene2/Gene3 region (Gene2 and Gene3 are overlaped).
• Gene4-Gene5-part1 : 1st block which have a common LD block overlapping both Gene4 and Gene5 gene regions.
• Gene4-Gene5-part2 : 2nd block which have a common LD block overlapping both Gene4 and Gene5 gene regions.
• inter-Gene5:Gene6-part1 : 1st block located at the intergenic region between Gene5 and Gene6
• Gene6-part1 : 1st block overlapping Gene6 region
• after-Gene100-part1 : 1st block located at the region after the lastly located Gene100

2.4.3.3 Details

As with BigLD, GPART can run the algorithm based on LD measure, D’ by setting LD = "Dprime".

# gene based model
Gres1 = GPART(geno = geno, SNPinfo = SNPinfo, geneinfo = geneinfo)
# gene based model using LD measure Dprime
Gres2 = GPART(geno = geno, SNPinfo = SNPinfo, geneinfo = geneinfo, LD = "Dprime")

There are three GPART algorithms:

1. “Gene based”. first merge the overlapping LD blocks and gene regions, split the big LD blocks and merged gene-block regions, and then merge small consecutive blocks such as singletons to satisfy the pre-defined min-size threshold (set GPARTmode="geneBased", default).

2. “LDblock based only”. This algorithm does not use any gene information. It partitions the SNP sequence by splitting or merging the LD blocks to satisfy the given size threshold (GPARTmode="LDblockBased", and Blockbasedmode="onlyBlocks", default for LD block based methods).

3. “LDblock based - use gene info”. This use gene information only for the merging consecutive small blocks procedure. If consecutive small blocks of which size are less than threshold do overlap with gene region, we merge them as large as possible (at most max-size threshold), otherwise, we merge them as small as possible. (set GPARTmode="LDblockBased", and Blockbasedmode="useGeneRegions")

# LD block based - use only LD block information
Gres3 = GPART(geno = geno, SNPinfo = SNPinfo, geneinfo = geneinfo, 
              GPARTmode = "LDblockBased", Blockbasedmode = "onlyBlocks")
# LD block based - use gene information to merge singletons
Gres4 = GPART(geno = geno, SNPinfo = SNPinfo, geneinfo = geneinfo, 
              GPARTmode = "LDblockBased", Blockbasedmode = "useGeneRegions")

You can load gene information from a text file using the option geneinfofile.

# you can load text file containing gene information.
Gres5 = GPART(geno = geno, SNPinfo = SNPinfo, geneinfofile = "geneinfo.txt")

You can also load the gene information data from “Ensembl” or “UCSC” database via internal sub-function of GPART. You need to set the database and assembly via the options of function, geneDB for database and assembly for assembly. Default for geneDB is "ensembl" and default for assembly is "GRCh38".

When you use “ensembl” data, you can set the version of the data by using the option ensbversion. By default, we use “hgnc symbol” for the gene name. You can change it by the option geneid depending on the gene DB you choose. If there are more than 2 gene information with the same gene name, we union the regions with the same name.

# use `geneDB` option instead of `geneinfo` or `geneinfofile`
# default geneDB is "ensembl" and  default assembly is "GRCh38"
Gres6 = GPART(geno = geno[,1:1000], SNPinfo = SNPinfo[1:1000,],
              BigLDresult = res1, geneDB = "ensembl", assembly = "GRCh37")
Gres7 = GPART(geno = geno[,1:1000], SNPinfo = SNPinfo[1:1000,], 
              BigLDresult = res1, geneDB = "ucsc", assembly = "GRCh37" )