2011-11-25 45 views
0

我的数据集dart是一个尺寸为1981 x 278的矩阵。第一列包含从1到21的染色体数字,第二列是标记名称,第三列是距离(CM )。循环相同的步骤来创建一个图

下面的代码绘制了一条染色体的LD衰减。我想对21条染色体重复同样的事情(循环它们)。

任何帮助或评论将不胜感激。

dart<- read.csv("dartnonaR.csv") 
chr1 <- which(dart[, 1] == 1); 
mpos <- dart[chr1,2:3 ]; 
head(mpos); 
dart1 <- dart[chr1,]; 
dim(dart1); 
dart2 <- dart1[,-c(1,2,3)]; 
dart2 <- t(dart2); 
r2 <- (cor(dart2))^2; 
rownames(r2) <- mpos$MARKERS; 
mark <- rownames(r2); 
r2a <- r2; 
r2v <- NULL; 
distance <- NULL; 

for(i in 1:144){ 
    for (j in (i+1):145){ 
     r2v <- c(r2v, r2a[i,j]) 
     distance <- c(distance, abs(mpos[mpos$MARKERS == mark[i],2] - mpos[mpos$MARKERS == mark[j],2])) 
     cat(i,j,"\n") 
    } 
}; 
plot(distance, r2v, xlab = "Distance in cM", ylab = "LD in r2"); 
+0

你能发表'str(dart)'的结果吗?或者甚至更好,组成一个模拟数据集(请参阅http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)。 –

回答

2

您需要在生成chr1子集的那一刻开始染色体上的循环。

要循环所有的染色体,你可以试试这个。我改编了一下你的代码。

dart <- read.csv("dartnonaR.csv") ## read data 
    savepdf = TRUE 
    for (k in 1:21){ ## start loop over chromosomes 
    chr <- which(dart[, 1] == k); ## assign data from col 1 to chr if equal to k 
    mpos <- dart[chr, 2:3 ]; ## create mpos 
    dart_chr <- dart[chr, ]; ## create dart_chr from dart 
    dart_chr2 <- t(dart_chr[, -c(1, 2, 3)]); ## get genomic data and transpose 
    r2 <- (cor(dart_chr2))^2; ## calculate r-square data 
    rownames(r2) <- mpos$MARKERS; ## Add rownames based on marker names 
    r2v <- NULL; ## initialize values 
    distance <- NULL; ## initialize values 
    for(i in 1:length(r2[,1])){ 
     for (j in (i+1):length(r2[1,])){ ## probably, could also be length(r2[1,]) + 1 , I'm not sure. 
     r2v <- c(r2v, r2[i, j]) 
     distance <- c(distance, abs(mpos[mpos$MARKERS == rownames(r2)[i], 2] - mpos[mpos$MARKERS == rownames(r2)[j], 2])) 
     cat(i, j, "\n") 
     } 
    }; 
    if(savepdf){ 
     pdf(file = paste('ld_decay_chr',k,'.pdf', sep = '')) 
     plot(distance, r2v, xlab = "Distance in cM", ylab = "LD in r2", main = paste('LD Decay chromosome', k)); 
     dev.off() 
    } 
    if(!savepdf){ 
     plot(distance, r2v, xlab = "Distance in cM", ylab = "LD in r2", main = paste('LD Decay chromosome', k)); 
    } 
    } 
0

要将所有情节组合成一个情节,一定要看看ggplot。更具体地说,facet_wrap和facet_grid函数。这些允许每个数据类别制作相同的图表,并将它们排列在一个图表格中。结合这些图可以轻松比较各类别之间的现货趋势。

查看http://had.co.nz/ggplot2/facet_grid.html的例子,包括漂亮的图片:)。

+0

与使用plot()相比,这还需要相当少的自定义代码。 –