2011-11-25 48 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()相比,這還需要相當少的自定義代碼。 –