0

我有一个不同染色体上的数字数据框（遗传学数据），可以将其视为分隔数字的因素。它看起来像这样（相邻的列包含每个位置的样本信息）：

awk '{print $2 "\t"   $3}' log_values | head 
Chr   Start    Sample1     Sample2
1   102447376  0.46957632  0.38415043
1   102447536  0.49194950  0.30094824
1   102447366  0.49874880 -0.17675325
2   102447366 -0.01910729  0.20264680
1   108332063 -0.03295081  0.07738970
1   109472445  0.02216355 -0.02495788

我想要做的是制作一系列从该文件中的其他列获取值的图。如果开始列中的值彼此足够接近，我想绘制覆盖范围的图，而不是为每一行绘制一个（这将代表不同区域和/或不同样本中的结果）。首先，如果 Start 列中的三个值彼此相距 1000 以内，我想绘制一个图。也就是说，从 A 到 B 到 C 包含 1000，因此 A 到 B <= 1000 和 B 到 C <= 1000 但 A 到 C 不必 <= 1000。在下面的代码中，这个 1000 是“ CNV_size”。“flanking_size”变量只是将绘图缩小一点，以便我可以给它一些上下文。

取样本值，第 1 行、第 2 行和第 3 行将突出显示为 Sample1 的一个图。这些样本数是 log2Ratios，所以我只想绘制重要的数。我将其定义为高于 0.4 或低于 -0.6。这意味着相同的三行不会产生样本 2 的图。

由于 Chr 列号/因子不同，因此不包括第四行。这是每列的单独图，仅显示满足此条件的行中的值。所以我可以为每个样本绘制多个图，但每组符合此标准的区域都将绘制在所有样本中。如果这没有意义，也许我下面的无效尝试将有助于解释我在胡说八道。

pdf("All_CNVs_closeup.pdf")

CNV_size <- 1000 # bp 
flanking_size <- 1000 # bp

#for(chr in 1:24){
for(chr in 1:1){
 #for(array in 1:24) {
for(array in 1:4) {

 dat <- subset(file, file$Chr == chr ) 
 dat <- subset(dat, dat[,array+6] > 0.4 | dat[,array+6] < -0.6) 
 if(length(dat$Start) > 1 ) { 
 dat <- dat[with(dat, order(Start)), ]

x=dat$Start[2:length(dat$Start)]-dat$Start[1:(length(dat$Start)-1)]
cnv <- 1
while(cnv <= length(x)) { 
for(i in cnv:length(x) ) {
    if(x[i] >= CNV_size) {
        plot_title <- paste(sample_info$Sample.ID[array], files[array],   sep = "   ") 
        plot(dat$Start, -dat[,array+6], main = plot_title , ylim = c(-2,2),      xlim = c(dat$Start[cnv] - flanking_size , dat$Start[i ] + flanking_size) , xlab = chr, ylab = "Log2 Ratio") 
        abline(h = 0.4, col="blue") 
        abline(h = 0, col="red") 
        abline(h = -0.6, col="blue") 
    break
    } # if(x[i] >= CNV_size) {      
    #if(x[i] < CNV_size) i <- i + 1
}   # for(i in cnv:length(x) ) {

cnv <- i 

} # while(x[cnv] <= length(x)) { 
  } # if(length(dat$Start) > 1 ) { 
  } # for(array in 1:24) {
} # for(chr in 1:24){



 dev.off()
4

1 回答 1

0

您可以编写一个循环来累积给定条件的索引，然后绘制每个窗口：

# Assuming your dataframe is called SNPs.
getChrWindows <- function(snps, windowSize) {
  curWindow <- 1
  windows <- list()
  for(i in 2:nrow(snps)) {
    # If pair is on the same chromosome and within the window, add it
    if(snps$Chr[i-1] == snps$Chr[i] && 
       (snps$Start[i] - snps$Start[i-1]) <= windowSize) {
      tryCatch({
        windows[[curWindow]]
      }, error = function(e) {
        windows[[curWindow]] <- c(i-1)    
      }, finally {
        windows[[curWindow]] <- c(windows[[curWindow]], i)
      })
    } else {
      # If there is an existing window, create a new one.
      tryCatch({
        windows[[curWindow]]
      }, error = function(e) {
        curWindow <- curWindow + 1
      })
    }
  }
  return(windows)
}

现在您可以获取 中所有窗口的列表data.frame，并绘制每个窗口：

windows <- getChrWindows(snps, 1000)
for (i in seq_along(windows)) {
  # plot snps[windows[[i]],] using your plotting code.
}
于 2013-10-07T05:26:44.827 回答