0

My data set dart is a matrix with dimensions 1981 x 278. The first column contains chromosome numbers from 1 to 21, the second column is the marker names and the third is the distances (CM).

The code below plots the LD decay for one chromosome. I want to repeat the same thing for the 21 chromosomes (looping over them).

Any help or comment will be appreciated.

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");
joran
  • 169,992
  • 32
  • 429
  • 468
hema
  • 725
  • 1
  • 8
  • 20
  • Can you post results of `str(dart)`? Or even better, make up a mock data set (see http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – Roman Luštrik Nov 25 '11 at 07:24

2 Answers2

2

You need to begin the loop over the chromosomes at the moment you generate the subset for chr1.

To loop over all the chromosomes, you could try this. I adapted you code a bit.

  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)); 
    }
  }
Mischa Vreeburg
  • 1,576
  • 1
  • 13
  • 18
0

To combine all plots into one plot, be sure to take a look at ggplot. More specifically the facet_wrap and facet_grid functions. These allow one to make identical plots per category of data, arranging them in a lattice of plots. Combining the plots allows easy comparison and spot trends between the categories.

See http://had.co.nz/ggplot2/facet_grid.html for examples, including nice pictures :).

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149