2楼: Originally posted by FMStation at 2019-02-12 02:08:49
##读入gff文件
header.lines <- readlines("plasmodb-33_pvivaxp01.gff", n = 30)

## 对数据清洗，获取该基因组每条染色体的长度大小，类似于产生一个染色体位置文件
ll <- h ...

##读入gff文件
header.lines <- readlines("plasmodb-33_pvivaxp01.gff", n = 30)

## 对数据清洗，获取该基因组每条染色体的长度大小，类似于产生一个染色体位置文件
ll <- header.lines[grepl(header.lines, pattern = "##sequence-region pvp01")]
gg <- data.frame(do.call(rbind, strsplit(ll, split = " ")))
gg[,3] <- as.numeric(as.character(gg[,3]))
gg[,4] <- as.numeric(as.character(gg[,4]))

##这里利用一个函数，将数据框转化成granges 对象
pvp01.genome <- togranges(gg[,c(2,3,4)])

##描绘基因组核型
kp <- plotkaryotype(genome=pvp01.genome)
kp <- plotkaryotype(genome=pvp01.genome, ideogram.plotter = null)
kpaddcytobandsasline(kp)

##读取该gff文件
features <- import("plasmodb-33_pvivaxp01.gff")

##将含有基因的行提取出来
genes <- features[features$type=="gene"]

##根据基因位置将基因描绘在染色体上
kpplotregions(kp, data=genes,col="blue")