OpenMendel系列-ADMIXTURE和OpenADMIXTURE

doc, genetics, julia, openmendel    
  1. ADMIXTURE.jl
  2. OpenADMIXTURE.jl
    1. 安装
    2. 用法
    3. 多线程
    4. GPU支持
    5. AdmixData

ADMIXTURE.jl

ADMIXTURE.jl ADMIXTURE.jl是ADMIXTURE的julia包装器, 等于直接运行ADMIXTURE再加了个Julia的壳, 就不赘述了, 直接学习OpenADMIXTURE.jl。

OpenADMIXTURE.jl

OpenADMIXTURE.jl

  • 是ADMIXTURE的julia重新实现, 采用更高效的多线程方案, 比原软件快8倍; 支持CUDA GPU计算; 输入为PLINK BED格式, 内存使用更小(比Float32,Float64小16x-32x)。

  • 假定个体不相关, 使用ML估计SNP的基因型数据祖先状态。支持利用sparse K-means via feature ranking (SKFR)方法选择SNP子集运行。

安装

  • Julia v1.6+

julia

# 安装软件包
using Pkg
pkg"add https://github.com/kose-y/SKFR.jl"
pkg"add https://github.com/OpenMendel/OpenADMIXTURE.jl"
# 安装示例应用相关包
pkg"add SnpArrays CSV DelimitedFiles"
# GPU支持
pkg"add CUDA"

julia

用法

julia

using LinearAlgebra, Random, SnpArrays
using OpenADMIXTURE
using CSV, DelimitedFiles

# 获取测试数据
filename = SnpArrays.datadir("EUR_subset.bed");

# 核心函数 run_admixture()
# 不使用SKFR:
d, clusters, aims = OpenADMIXTURE.run_admixture(filename, 4; T=Float64, use_gpu=false, rng=MersenneTwister(7856))
# 使用SKFR:
d, clusters, aims = OpenADMIXTURE.run_admixture(filename, 4; T=Float64, use_gpu=false, rng=MersenneTwister(7856), sparsity=10000, admix_rtol=1e-7, prefix="./EUR_subset")

julia
Note run_admixture()参数说明:

输入

  1. T: 估计精度, Float64或Float32, 默认是Float64;

  2. use_gpu: 是否使用GPU, 默认false;

  3. rng: 随机数生成方法, 默认Random.GLOBAL_RNG, 采用MersenneTwister是SKFR.jl包中推荐用法, 好像是可以在多线程中稳定随机;

  4. prefix: SKFR输出的SNP bed文件的前缀, 输出文件格式$(prefix)_$(k)_$(sparsity)aims.bed;

  5. sparsity: SKFR选择的SNP数量, 默认nothing(跳过SKFR);

  6. skfr_tries: 运行SKFR多次并选择最佳集群, 默认1;

  7. skfr_max_inner_iter: 运行每个SKFR最多迭代次数或直到收敛, 默认50;

  8. admix_n_iter: ADMIXTURE的最大迭代次数, 默认1000;

  9. admix_rtol: 根据对数似然的相对变化收敛标准, 默认1e-7;

  10. admix_n_em_iter: EM迭代次数, 用于获得最佳初始值, 默认5;

  11. Q: quasi-Newton acceleration中的步数(不知道干啥的), 默认3;

输出: (d,clusters,aims)

  • d: Admixture结果数据结构AdmixData, 有很多信息, 其中d.p存储AF, d.q存储admixture proportions;

  • clusters: 样本聚类标签, 如果sparsity == nothing则不输出;

  • aims: SKFR所选的SNP索引, 按照重要性降序排列, 如果sparsity == nothing则不输出;

多线程

默认使用julia占用的线程数, 所以启动julia的时候设置多线程: julia -t 8

GPU支持

use_gpu = true开启GPU, 将会把gradients and Hessians of the loglikelihood计算放到GPU中执行。

AdmixData

julia

AdmixData{T}(I, J, K, Q; skipmissing=true, rng=Random.GLOBAL_RNG)
Constructor for Admixture information.

Arguments:
- I: Number of samples
- J: Number of variants
- K: Number of clusters
- Q: Number of steps used for quasi-Newton update
- skipmissing: skip computation of loglikelihood for missing values. Should be kept `true` in most cases
- rng: Random number generation.
- 还有其他一堆...

julia