R语言Julia以及全基因组选择
最近在学Julia语言，想测试一下和R的区别，发现前辈的博客，翻译时不禁感慨，这是2018年了，博客是2010年的，8年已过，我才听说Julia。但……不晚！
文章来源： 祝福病人早日康复短语/r-julia-and-genome-wide-lection/
有一些琐事的事情，以及一些代码片段，这些是我在2010年参加了基因组选择的summer school（/wiki/pages/E4o0S0C7/Cour_materials.html），数据共有2000个个体，20000个SNP（0,1,2），然后使用MCMC 计算育种值，使用的是R语言。
```
nmarkers = 2000; # number of markersstartMarker = 1981; # t to 1 to u 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 u 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，特别的漂亮鸡今年多大
