我当时是这么解决的,就是把每一个重复区域的起始和终止位点都记录一下,然后画矩形。
效果图类似:
[attach]276[/attach]
代码类似:
[code]hg19_g_band = read.table("~/menghw_HD/reference/hg19_G-band.txt",header = F,sep = "\t")
plot.chromosome <- function(chrom_index,chrom_start=0,chrom_end=250e6,track_height=8,border.lwd = 1,df.peak = hg19_g_band ){
  # initialize
  # chrom_index = 1
  # chrom_start = 0
  # chrom_end = chrom.length(chrom_index)
  # 将chrom_index生成bed文件中用的chrom_name
  if(chrom_index <=22){
    chrom_name = paste0("chr",chrom_index)
  }else if(chrom_index == 23){
    chrom_name = "chrX"
  }else if(chrom_index == 24){
    chrom_name = "chrY"
  }
  
  # 只画某1条染色体,筛选数据
  bed_table = df.peak[df.peak[,1]==chrom_name,]
  filter_vector = (bed_table[,2] >= chrom_start & bed_table[,2] <= chrom_end) & (bed_table[,3] >= chrom_start & bed_table[,3] <= chrom_end)
  bed_table = bed_table[filter_vector,]
  
  # 只使用有用的数据信息
  track_start = bed_table[,2]
  track_end = bed_table[,3]
  track_value = bed_table[,5]
  track_info = bed_table[,6]
  ## use 255 gray pattern
  alpha_vector = floor(track_value / 600 * 255) # 这里设置600纯是为了好看
  alpha_vector[alpha_vector>255] = 255
  color_vector = rgb(t(col2rgb("black")),alpha = alpha_vector,maxColorValue = 255)
    
  # rect 坐标
  track_x1 <- as.vector(track_start)
  track_x2 <- as.vector(track_end)
  track_y1 <- rep((-1)*(track_height %/% 2),nrow(bed_table))
  track_y2 <- rep(track_height %/% 2,nrow(bed_table))
  
  track_y1[track_info=="acen"] <- (-1)*(track_height %/% 2)/2
  track_y2[track_info=="acen"] <- (track_height %/% 2)/2
  
  plot(x=c(chrom_start,chrom_end),y=c((-1)*(track_height %/% 2),(track_height %/% 2)),type="n",frame.plot = F,cex.axis=2,cex.lab=2,ylab = "",xlab = "",xaxt = "n",yaxt="n")
  rect(track_x1,track_y1,track_x2,track_y2,col = color_vector,border = T,lwd = border.lwd)
}[/code]然后hg19-G-band.txt文件的格式类似于:
[code]chr1  0 2300000 p36.33  0 gneg
chr1  2300000 5400000 p36.32  250 gpos25
chr1  5400000 7200000 p36.31  0 gneg
chr1  7200000 9200000 p36.23  250 gpos25
chr1  9200000 12700000  p36.22  0 gneg
chr1  12700000  16200000  p36.21  500 gpos50
chr1  16200000  20400000  p36.13  0 gneg
chr1  20400000  23900000  p36.12  250 gpos25
chr1  23900000  28000000  p36.11  0 gneg
chr1  28000000  30200000  p35.3 250 gpos25
chr1  30200000  32400000  p35.2 0 gneg
chr1  32400000  34600000  p35.1 250 gpos25
chr1  34600000  40100000  p34.3 0 gneg
chr1  40100000  44100000  p34.2 250 gpos25
chr1  44100000  46800000  p34.1 0 gneg
chr1  46800000  50700000  p33 750 gpos75
chr1  50700000  56100000  p32.3 0 gneg
chr1  56100000  59000000  p32.2 500 gpos50
chr1  59000000  61300000  p32.1 0 gneg
chr1  61300000  68900000  p31.3 500 gpos50[/code]
 
 
                                 
                                                                
                                     阅读全文
                                
                                
                                     收起全文