цикл по тем же шагам, чтобы создать сюжет - PullRequest
0 голосов
/ 25 ноября 2011

Мой набор данных dart - это матрица с размерами 1981 x 278. Первый столбец содержит номера хромосом от 1 до 21, второй столбец - это имена маркеров, а третий - расстояния (см).

Приведенный ниже код показывает распад 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");

Ответы [ 2 ]

2 голосов
/ 25 ноября 2011

Вам нужно начать цикл над хромосомами в тот момент, когда вы генерируете подмножество для 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 голосов
/ 27 ноября 2011

Чтобы объединить все графики в один график, обязательно посмотрите на ggplot. Более конкретно, функции facet_wrap и facet_grid. Это позволяет составлять идентичные графики для каждой категории данных, размещая их в решетке графиков. Объединение графиков позволяет легко сравнивать и определять тренды между категориями.

См. http://had.co.nz/ggplot2/facet_grid.html для примеров, включая красивые картинки:).

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...