Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

comment out QTLPlot module and update project.toml #43

Merged
merged 12 commits into from
Feb 13, 2024
Merged
8 changes: 2 additions & 6 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -10,23 +10,19 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LossFunctions = "30fc2ffe-d236-52d8-8643-a9d8f7c094a7"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"

[compat]
Conda = "1"
DelimitedFiles = "1"
Distributions = "^0.23, 0.24, 0.25"
LossFunctions = "~0.11"
PyPlot = "~2.11"
Revise = "3.5"
StaticArrays = "1.7"
StaticArrays = "1.9"
StatsBase = "0.33, 0.34"
Statistics = "1"
StatsBase = "0.33"
julia = "~1.9"

[extras]
Expand Down
7 changes: 1 addition & 6 deletions docs/src/functions.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,6 @@ Modules = [FlxQTL, flxMLMM, EcmNestrv]
Modules = [GRM]
```

## Visualization (QTLplot)

```@autodocs
Modules = [QTLplot]
```

## Multivariate Linear Models (MLM)

```@autodocs
Expand All @@ -35,3 +29,4 @@ Modules = [MLM]
```@autodocs
Modules = [Util]
```

6 changes: 4 additions & 2 deletions docs/src/guide/analysis.md
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,9 @@ third marker, type `B[:,:,3]`.

## Generating plots

To produce a plot (or plots) for LOD scores or effects, you need first a struct of arrays, `layers` consisting of chromosomes, marker positions,
The `QTLplot` module is currently unavailable but will replaced with [BigRiverQTLPlots.jl](https://github.com/senresearch/BigRiverQTLPlots.jl) soon.

<!-- To produce a plot (or plots) for LOD scores or effects, you need first a struct of arrays, `layers` consisting of chromosomes, marker positions,
LOD scores (or effects), which should be `Array{Float64,2}`. You can then generate one genome scan result or multiple genenome scan results on one plot. Note that the color is randomly selected to generate a plot.
The function `plot1d` has more keyword argument options: `yint=[]` for a vector of y-intercept(s), `yint_color=["red"]` for a vector of y-intercept
color(s), `Legend=[]` for multiple graphs, `loc="upper right"` for the location of `Legend`, etc.
Expand All @@ -155,7 +157,7 @@ plot1d(Arab_lod;title= "LOD for Arabidopsis thaliana : Fitness (2 site by 3 year
![arabidopsis](arab-lod.png)



-->

## Performing a permutation test

Expand Down
2 changes: 1 addition & 1 deletion docs/src/guide/tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Or, equivalently,
```julia
julia> using Pkg; Pkg.add("FlxQTL")
```
Currently Julia `1.5` supports for the package.
<!-- Currently Julia `1.5` supports for the package. -->


To remove the package from the Julia REPL,
Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@ structured multivariate traits.

- Genome scan (1D, 2D) for univariate, multivariate trait(s), and genotype (probability) data
- LOCO (Leave One Chromosome Out) support for genome scan
- Visualization (1D, 2D plots)
- Computation for Genetic (or Climatic) Relatedness matrix (or kinship)
- CPU parallelization
<!-- - Visualization (1D, 2D plots) -->

## Guide

Expand Down
16 changes: 8 additions & 8 deletions src/geneScan1D.jl
Original file line number Diff line number Diff line change
Expand Up @@ -70,36 +70,36 @@
end

## upgraded version:
function marker1Scan(m,q,kmin,cross,λg,Y1,Xnul_t,X1,Z1,Nullpar::Result,Σ1,τ2_0;ρ=0.001,tol0=1e-3,tol1=1e-4,nchr=0)

Check warning on line 73 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L73

Added line #L73 was not covered by tests

nmar=size(X1,1);
if (cross==1) ## scanning genotypes
B0=hcat(Nullpar.B,zeros(Float64,q));

Check warning on line 77 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L75-L77

Added lines #L75 - L77 were not covered by tests

lod=@distributed (vcat) for j=1:nmar
XX=vcat(Xnul_t,@view X1[[j],:])
B0,τ2,Σ,loglik0 =ecmLMM(Y1,XX,Z1,B0,τ2_0,Σ1,λg,ones(m);tol=tol0)
lod0= (loglik0-Nullpar.loglik)/log(10)
est1=ecmNestrvAG(lod0,kmin,Y1,XX,Z1,B0,τ2,Σ,λg,ones(m);ρ=ρ,tol=tol1,numChr=nchr,nuMarker=j)
[(est1.loglik-Nullpar.loglik)/log(10) est1]
end

Check warning on line 85 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L79-L85

Added lines #L79 - L85 were not covered by tests

else # cross>1
## scanning genotype probabilities

B0=hcat(Nullpar.B,zeros(Float64,q,cross-1))

Check warning on line 90 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L90

Added line #L90 was not covered by tests

lod=@distributed (vcat) for j=1:nmar
XX= vcat(Xnul_t, @view X1[j,2:end,:])
B0,τ2,Σ,loglik0 =ecmLMM(Y1,XX,Z1,B0,τ2_0,Σ1,λg,ones(m);tol=tol0)
lod0= (loglik0-Nullpar.loglik)/log(10)
est1=ecmNestrvAG(lod0,kmin,Y1,XX,Z1,B0,τ2,Σ,λg,ones(m);ρ=ρ,tol=tol1)
[(est1.loglik-Nullpar.loglik)/log(10) est1]
end

Check warning on line 98 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L92-L98

Added lines #L92 - L98 were not covered by tests

end

return lod[:,1],lod[:,2]

Check warning on line 102 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L102

Added line #L102 was not covered by tests
end


Expand Down Expand Up @@ -221,11 +221,11 @@
for i=1:nChr
maridx=findall(XX.chr.==Chr[i])
# Xnul_t=Xnul*Tg[:,:,i]';
@fastmath @inbounds Xnul_t=BLAS.gemm('N','T',Xnul,@view Tg[:,:,i])
@fastmath @inbounds Xnul_t=BLAS.gemm('N','T',Xnul, Tg[:,:,i])
if (cross!=1)
@fastmath @inbounds @views Y2,X1=transForm(Tg[:,:,i],Y1,X0[maridx,:,:],cross)
@fastmath @inbounds Y2,X1=transForm(Tg[:,:,i],Y1,X0[maridx,:,:],cross)
else
@fastmath @inbounds @views Y2,X1=transForm(Tg[:,:,i],Y1,XX.X[maridx,:],cross)
@fastmath @inbounds Y2,X1=transForm(Tg[:,:,i],Y1,XX.X[maridx,:],cross)
end
#parameter estimation under the null
est00=nulScan(init,1,Λg[:,i],λc,Y2,Xnul_t,Z1,Σ1;ρ=ρ,itol=itol,tol=tol)
Expand Down Expand Up @@ -297,11 +297,11 @@
for i=1:nChr
maridx=findall(XX.chr.==Chr[i])
# Xnul_t=Xnul*Tg[:,:,i]';
@fastmath @inbounds Xnul_t=BLAS.gemm('N','T',Xnul,@view Tg[:,:,i])
@fastmath @inbounds Xnul_t=BLAS.gemm('N','T',Xnul,Tg[:,:,i])
if (cross!=1)
@fastmath @inbounds @views Y2,X1=transForm(Tg[:,:,i],Y1,X0[maridx,:,:],cross)
@fastmath @inbounds Y2,X1=transForm(Tg[:,:,i],Y1,X0[maridx,:,:],cross)
else
@fastmath @inbounds @views Y2,X1=transForm(Tg[:,:,i],Y1,XX.X[maridx,:],cross)
@fastmath @inbounds Y2,X1=transForm(Tg[:,:,i],Y1,XX.X[maridx,:],cross)
end
#parameter estimation under the null
est00=nulScan(init,1,Λg[:,i],λc,Y2,Xnul_t,Σ1;ρ=ρ,itol=itol,tol=tol)
Expand Down Expand Up @@ -363,9 +363,9 @@
# Xnul_t=Xnul*Tg[:,:,i]';
@fastmath @inbounds Xnul_t=BLAS.gemm('N','T',Xnul,@view Tg[:,:,i])
if (cross!=1)
@fastmath @inbounds @views Y,X=transForm(Tg[:,:,i],Y0,X0[maridx,:,:],cross)
@fastmath @inbounds Y,X=transForm(Tg[:,:,i],Y0,X0[maridx,:,:],cross)
else
@fastmath @inbounds @views Y,X=transForm(Tg[:,:,i],Y0,XX.X[maridx,:],cross)
@fastmath @inbounds Y,X=transForm(Tg[:,:,i],Y0,XX.X[maridx,:],cross)
end
#parameter estimation under the null
est00=nulScan(init,1,Λg[:,i],Y,Xnul_t;itol=itol,tol=tol,ρ=ρ)
Expand Down Expand Up @@ -414,14 +414,14 @@
# export geneScan,marker1Scan
## upgraded version: H0:MLMM, H1:FlxQTL

function obtainKc(T::Array{Float64,2},λ::Array{Float64,1},Y0::Array{Float64,2},Z0::Array{Float64,2};

Check warning on line 417 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L417

Added line #L417 was not covered by tests
Xnul::Array{Float64,2}=ones(1,size(Y0,2)),itol=1e-2,tol::Float64=1e-3,ρ=0.001)

# T,λ = K2eig(Kg)
init=initial(Xnul,Y0,Z0,false)
Y, Xnul_t=transForm(T,Y0,Xnul,1)
est0=nulScan1(init,1,λ,Y,Xnul_t,Z0;itol=itol,tol=tol,ρ=ρ)
return est0.Vc, Y, Xnul_t

Check warning on line 424 in src/geneScan1D.jl

View check run for this annotation

Codecov / codecov/patch

src/geneScan1D.jl#L421-L424

Added lines #L421 - L424 were not covered by tests
end

include("geneScan1.jl")
Expand Down
2 changes: 1 addition & 1 deletion src/transformation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@
else #cross>1 size(X0)= (p,cross,n)
p=size(X0,1); n=size(X0,3)
X=zeros(p,cross,n)
@fastmath @inbounds @views for j=1:cross
@fastmath @inbounds @views for j=1:cross
# X[:,j,:] = X0[:,j,:]*Tg'
X[:,j,:]= BLAS.gemm('N','T',X0[:,j,:],Tg)
end
Expand Down Expand Up @@ -208,8 +208,8 @@
return Init(init_val.B',τ2,init_val.Σ)
else #H0 for MLMM

lmul!(sqrt(1/m),init_val.Σ)
return Init0(init_val.B',0.5*init_val.Σ,0.5*init_val.Σ)

Check warning on line 212 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L211-L212

Added lines #L211 - L212 were not covered by tests
end

end
Expand Down Expand Up @@ -242,10 +242,10 @@
return nulpar
end

function nulScan(init::Init1,kmin,λg,λc,Y1,Xnul_t,Z1;ρ=0.001,itol=1e-3,tol=1e-4)

Check warning on line 245 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L245

Added line #L245 was not covered by tests

B0,τ2_0,Σ1,loglik0 =ecmLMM(Y1,Xnul_t,Z1,init.B,init.τ2,init.Σ1,λg,λc;tol=itol)
nulpar=NestrvAG(kmin,Y1,Xnul_t,Z1,B0,τ2_0,Σ1,λg,λc;ρ=ρ,tol=tol)

Check warning on line 248 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L247-L248

Added lines #L247 - L248 were not covered by tests

return nulpar

Expand All @@ -261,16 +261,16 @@
end

#new version to estimate Kc
function nulScan1(init::Union{Init1,Init0},kmin,λg,Y1,Xnul_t,Z;ρ=0.001,itol=1e-3,tol=1e-4)

Check warning on line 264 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L264

Added line #L264 was not covered by tests

if (typeof(init)==Init1)
B0,Kc_0,Σ1,loglik0 =ecmLMM(Y1,Xnul_t,Z,init.B,init.Kc,init.Σ1,λg;tol=itol)
nulpar=NestrvAG(kmin,Y1,Xnul_t,Z,B0,Kc_0,Σ1,λg;tol=tol,ρ=ρ)

Check warning on line 268 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L266-L268

Added lines #L266 - L268 were not covered by tests
else #typeof(Init)==Init0)
B0,Kc_0,Σ1,loglik0 =ecmLMM(Y1,Xnul_t,Z,init.B,init.Vc,init.Σ,λg;tol=itol)
nulpar=NestrvAG(kmin,Y1,Xnul_t,Z,B0,Kc_0,Σ1,λg;tol=tol,ρ=ρ)

Check warning on line 271 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L270-L271

Added lines #L270 - L271 were not covered by tests
end
return nulpar

Check warning on line 273 in src/transformation.jl

View check run for this annotation

Codecov / codecov/patch

src/transformation.jl#L273

Added line #L273 was not covered by tests
end


Expand Down
Loading