我当时是这么解决的,就是把每一个重复区域的起始和终止位点都记录一下,然后画矩形。
效果图类似:
[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]
阅读全文
收起全文