全基因组分析做图函数

很不错的,值得拥有

不同的语言，有不同的风格。R是有贵族气息的，她是一个指挥着，而不是执行者。就像作战指挥的将军，你让他去冲锋陷阵，当然也不是不行，但是毕竟年龄比较老了，行动不便。与年轻的士兵比较，自然逊色不少了。大将自然要发挥大将的优点，指挥需要的是那种智慧和狡黠。而不是极强的执行力。需要执行某个任务的时候，让身边的士兵（C,Fortran,C++等）去就好了。

现在举一个例子来说明。

假设有100个教室，统计每个教室的人数。士兵（C）来做的话，跑到第一教室，数一数第一个教室的人数。跑到第二个教室去，数一数第二个教室的人数。。。。。跑第100个教室，数一数第100个教室的人数。好在年轻士兵执行力强，这点体力还是有的，不用休歇。让将军去做（R),如果你也说跟士兵一样，一个一个跑过去，老胳膊老腿了，自然体力不行，干一阵得歇一阵，顶不住啊。所以就老有人抱怨，唉，这个老不死的有啥用啊，赶紧把它给换了……。但是如果你说，让R叫100个士兵分别去100个教室统计人数，然后回来报告给R，统计总人数，这下不用跑腿了，而且还是齐上阵啊，那效率，刷刷的。。。。所以也有人说，这么NB的将军你就还说要换，太不要脸了。。。

当然，也有的时候，非得要一步一步地执行才可以。比如砌墙，你不能说我从上面砌，你从下面砌，另外一个人从中间砌。这个时候英雄也有悲哀的时候。咋办啊？这个时候啊，就想个办法，先把砖块搞大，搞成一大块一大块的，然后再来砌，三下五除二，比你蹲在那里，沟着要，辛苦地一小块一小块磨着要爽多了吧。

所以，如果用R的人，千万别把他当成士兵，而要把他想成将军，他是指挥者，需要智慧去调动他的千军万马，而不是让他一个人去冲锋陷阵。现在来举一个小例子来说明那高贵的血统，是何等的优雅～。

假设我有一个数组如下：

很不错的,值得拥有

一共有7组因素，每个因素内有4个水平。你要求出每个水平内和的平方，然后再把各个水平的结果求和。

比如，对第一组因素，第一个水平的结果是（0.112+0.09+0.47+0.38）^2,第二个水平上（0.67+0.03-0.58-0.34）^2,以此类推，求出所有的四个水平内和的平方，然后再求和。遍历这7个因素。 方法一：土士兵的方式：

A<-matrix(sample(1:4,100000,replace=T),100,1000)

phe<-rnorm(100)

score<-rep(0,1000)

for (i in 1:1000){

group<-rep(0,4)

for (j in 1:100){

group[A[j,i]]=group[A[j,i]]+phe[j]

}

score<-sum(group)

}

这个方式，估计是学C的最擅长的了。

方法二：稍微优雅一点的方式

for (i in 1:1000){

score<-lapply(list(A[,i]),function(x,phe){sum(sum(phe)^2))}

}

这个，是比较好一点点的方法

当然，你会用更优雅一点点的方法

很不错的,值得拥有

score<-apply(A,2,phe,function(x){lapply(list(x),sum(sum(phe)^2))} 一句话完成了，爽歪歪。不过你也会发现，其实它的执行速度不太爽。 这个时候狡黠的将军就有了狡黠的歪点子了，他说，不行，咱得多人并行作业，这样太慢了。 所以，你开始了如下的代码： group<-matrix(100,4) score<-rep(0,1000) for (i in 1:1000){ group[,A[,i]]<-phe score<-sum(rowSum(group)^2) } 这样，实现了R的并行计算，它的坐镇指挥功能得到发挥。当然，还有更优雅的绝招还在后头，它甚至优雅地胜出了以Fortran这类计算语言 的尖兵。 优雅而高贵的R，千万别当士兵使唤。 如何用R写一个计算等位基因频率的函数

如果有一个SNP芯片的数据，要计算每个标记的等位基因频率，如何用R写这么一个函数，使它的计算速度达到极限？

假定：1000个个体，60,000 markers

SNP基因型为1 1， 1 2，2 2 的形式。请大家自由发挥

~~~~~~~~~~~~~~~~~~~~~~~~

allele_freq<-function(genotype){

nSNP<-ncol(genotype)

nind<-nrow(genotype)

g1<-genotype[,seq(1,nSNP,2)]

g2<-genotype[,seq(2,nSNP,2)]

freq<-(colSums(g1-1)+colSums(g2-1))/nind/2

freq

}

#test

geno<-matrix(sample(1:2,1000*12000,replace=T),nrow=1000) #simulate 1000 inpiduals # and 60,000SNPs

myfreq<-allele_freq(geno)

~~~~~~~~~~~~~~~~~~~~~·

哈哈，长啸一声，再次更新算法，真是太佩服我自己了。。。。

以小马的例子来算

t<-matrix(sample(c(11,12,21,22),30,T),nrow=10,ncol=3)

# t<-matrix(sample((c(11,12,21,22),1.2e8,T),nrow=1e3,ncol=6e4)

a<-1:3

t<-as.data.frame(t)

names(t)<-paste("M",a,sep="_")

#new method to calculate allele freq

freq<-function(geno){

#left hand allele: genotype 11 return 0, 12 return 0, 21 return 1, 22 return 1

allele_1<-geno%/%10-1

#right hand allele: genotype 11 return 0, 22 return 0, 12 return 1, 22 return 1

allele_2<-geno%%10-1

rm(geno);gc();

很不错的,值得拥有

alle_2<-colMeans(allele_1+allele_2)/2 geno_22<-colMeans(allele_1*allele_2) geno_12<-colMeans((allele_1-allele_2)^2) myres<-cbind(alle_2,geno_12,geno_22) myres } freq(t)

#this Function used for draw the figure for GWAS analysis

# data consturcted by atleast 3 columns,P-value,chromosome,position # correspond names of variables are pval CHR position'

# the names of the variables are very important, you cound change the names #chrpos=T option means the map distance was measured by each chromosome, #which need to be transform to genomewide distance by adding each distance # together

# if each chromosome distance was starts from 0 cM, then the plot

##bugs fixed

#2010-12-10 update for automatic adding each chromosome together

#2011-1-20 update for not sorted chromosome and position,with options sort=T

############################################################## ### FUNCTION plot.GWAS() ###

##############################################################

"plot.GWAS" <- function(x, y, ..., ystart = 0, log = F, chrpos=F,

ylim, delta = 1,sort=T) {

x<-data.frame(x)

Pv = x$pval

if (sort==T)x<-x[do.call(order,x[,c('CHR','position')]),]

if (log == T) Pv = -log10(x$pval)

mxlog <- ceiling(max(Pv, na.rm = T))

if (!missing(ylim)) {

ylim <- ylim

}

else {

ylim <- c(ystart, mxlog)

}

# detected if the distances were chromosome specific

if ((leng …… 此处隐藏：3653字，全部文档内容请下载后查看。喜欢就下载吧 ……