最近在学Julia语言，想测试一下和R的区别，发现前辈的博客，翻译时不禁感慨，这是2018年了，博客是2010年的，8年已过，我才听说Julia。但……不晚！
文章来源： https://www.r-bloggers.com/r-julia-and-genome-wide-selection/
有一些琐事的事情，以及一些代码片段，这些是我在2010年参加了基因组选择的summer school（http://taurus.ansci.iastate.edu/wiki/pages/E4o0S0C7/Course_materials.html），数据共有2000个个体，20000个SNP（0,1,2），然后使用MCMC 计算育种值，使用的是R语言。
```
nmarkers = 2000; # number of markersstartMarker = 1981; # set to 1 to use allnumiter = 2000; # number of iterationsvara = 1.0/20.0;
# input datadata = matrix(scan("trainData.out0"),ncol=nmarkers+2,byrow=TRUE);
nrecords = dim(data)[1];
beg = Sys.time()# x has the mean followed by the markersx = cbind(1,data[,startMarker:nmarkers]);y = data[,nmarkers+1];
a = data[,nmarkers+2];# inital valuesnmarkers = nmarkers - startMarker + 1;
mean2pq = 0.5; # just an approximationscalea = 0.5*vara/(nmarkers*mean2pq); # 0.5 = (v-2)/v for v=4size = dim(x)[2];
b = array(0.0,size);
meanb = b;
b[1] = mean(y);
var = array(0.0,size);# adjust y
ycorr = y - x%*%b; # MCMC samplingfor (iter in 1:numiter){ # sample vare
vare = ( t(ycorr)%*%ycorr )/rchisq(1,nrecords + 3); # sample intercept
ycorr = ycorr + x[,1]*b[1];
rhs = sum(ycorr)/vare;
invLhs = 1.0/(nrecords/vare);
mean = rhs*invLhs;
b[1] = rnorm(1,mean,sqrt(invLhs));
ycorr = ycorr - x[,1]*b[1];
meanb[1] = meanb[1] + b[1]; # sample variance for each locus
for (locus in 2:size){
var[locus] = (scalea*4+b[locus]*b[locus])/rchisq(1,4.0+1)
}# sample effect for each locus
for (locus in 2:size){ # unadjust y for this locus
ycorr = ycorr + x[,locus]*b[locus];
rhs = t(x[,locus])%*%ycorr/vare;
lhs = t(x[,locus])%*%x[,locus]/vare + 1.0/var[locus];
invLhs = 1.0/lhs;
mean = invLhs*rhs;
b[locus]= rnorm(1,mean,sqrt(invLhs)); #adjust y for the new value of this locus
ycorr = ycorr - x[,locus]*b[locus];
meanb[locus] = meanb[locus] + b[locus];
}
}
Sys.time() - beg
meanb = meanb/numiter;
aHat = x %*% meanb;```
我使用Julia去做类似的事情：需要定义几个新的变量，将基因组数据，表型数据以及育种值数据读进矩阵里面，写几个循环，进行向量的运算。
```
nmarkers = 2000 # Number of markersstartmarker = 1981 # Set to 1 to use allnumiter = 2000 # Number of iterationsdata = dlmread("markers.csv", ',')
(nrecords, ncols) = size(data)
tic()#this is the mean and markers matrixX = hcat(ones(Float64, nrecords), data[:, startmarker:nmarkers])y = data[:, nmarkers + 1]a = data[:, nmarkers + 2]nmarkers = nmarkers - startmarker + 1vara = 1.0/nmarkersmean2pq = 0.5scalea = 0.5*vara/(nmarkers*mean2pq) # 0.5 = (v-2)/v for v=4ndesign = size(X, 2)b = zeros(Float64, ndesign)meanb = zeros(Float64, ndesign)
b[1] = mean(y)varian = zeros(Float64, ndesign)# adjust yycorr = y - X * b # MCMC samplingfor i = 1:numiter # sample vare
vare = dot(ycorr, ycorr )/randchi2(nrecords + 3) # sample intercept
ycorr = ycorr + X[:, 1] * b[1];
rhs = sum(ycorr)/vare;
invlhs = 1.0/(nrecords/vare);
mn = rhs*invlhs;
b[1] = randn() * sqrt(invlhs) + mn;
ycorr = ycorr - X[:, 1] * b[1];
meanb[1] = meanb[1] + b[1]; # sample variance for each locus
for locus = 2:ndesign
varian[locus] = (scalea*4 + b[locus]*b[locus])/randchi2(4.0 + 1);
end # sample effect for each locus
for locus = 2:ndesign # unadjust y for this locus
ycorr = ycorr + X[:, locus] * b[locus];
rhs = dot(X[:, locus], ycorr)/vare;
lhs = dot(X[:, locus], X[:, locus])/vare + 1.0/varian[locus];
invlhs = 1.0/lhs;
mn = invlhs * rhs;
b[locus] = randn() * sqrt(invlhs) + mn; #adjust y for the new value of this locus
ycorr = ycorr - X[:, locus] * b[locus];
meanb[locus] = meanb[locus] + b[locus];
end
end
toc()meanb = meanb/numiter;aHat = X * meanb;```
这两个代码比较相似，但是也有一些不同：
第一个读入的数据是二进制的，我不知道Julia如何操作，所以我就将其转为csv，然后读取。
为了防止名称重复，R中可以随意命名，但是Julia不行，所以我在Julia程序中进行了修改。
R中可以赋值，a=b，你改变a和b都没问题，但是Julia中你动了b，a也变了。语法不太一样。
Julia中向量和数组操作不太一样。
因为我的代码太粗糙了，Julia的代码在速度上不是很突出。当我变化了marker的数目（数据量更大），发现Julia的运算速度大约是R的2.8倍。Julia在官网上宣称，在进行数值运算时，速度是R的100倍，但是我的结果没有达到这么高。
在1996年或者1997年是，我由SAS转到了ASReml进行基因组数据分析，它大约提高了1~2倍的速度，而且支持了更多的模型。
现在，又到了更换软件的时候了，由R转到Julia，特别是基因组选择方面，Julia性能非常优秀。
看一下Julia的logo，特别的漂亮！！！
