没空改了，等有空了再来修改格式吧
############genomicranges的学习################

首先使用granges来创建一个对象

gr <- GRanges(
seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
score = 1:10,
GC = seq(1, 0, length=10))

gr

查看对象中的内容

seqnames(gr) #查看基因的名称
ranges(gr) #查看基因的区域
strand(gr) #查看基因位于那条链
start(gr)
end(gr)
width(gr)
names(gr)
length(gr)

任何文件系统中的数据分为数据和元数据

数据是指普通文件中的实际数据，

而元数据指用来描述一个文件的特征的系统数据，诸如访问权限、文件拥有者以及文件数据块的分布信息等等，

而在此处元数据就是指基因的附加信息

granges(gr) #查看基因基本信息
as.data.frame(granges(gr))#将基因的基本信息转换为数据框

mcols(gr) #查看基因的元信息
as.data.frame(mcols(gr))#将基因的元数据转换为数据框

如果要查看单独的元数据，可以使用$

grGC
mcols(gr)GC

切割和合并GRanges

sp<-split(gr,rep(1:2,each=5))
c(sp[[1]],sp[[2]])

GRanges的亚集合

gr[2:3]
gr[2:3,"GC"]#后面的是GRanges的元数据列

将某一个基因进行重复

singles<-split(gr,names(gr))#将gr根据names来分割
grMod<-gr[1:4]
grMod
grMod[2]<-singles[[1]]
grMod

grMod[2]<-grMod[1]

repeat, reverse, 重复和反向

rep(singles[[2]],times=3)
rev(gr)

select,截取

head(gr,n=2)
tail(gr,n=2)
window(gr,start=2,end=4)
gr[IRanges(start=c(2,7), end=c(3,9))]

基本的一些操作

g<-gr[1:3]
g<-append(g,singles[[10]])
g

对于GRanges的方法，分为内部，外部，之间

内部flank, resize, shift

flank(g,10)#取得上游的10bp
flank(g,10,start=F)#取得下游的10bp
shift(g,5)#向上游偏移
shift(g,-5)#向下游偏移
resize(g,30)#调整width，只调整下游

如果移动过有负数，可以将其变为1

flank(g,150)
start(g[start(g)<1])<-1

外部reduce，gaps，range，disjoin，coverage

reduce会对区间进行合并overlap，得到一个简化的区间

reduce(g)
gaps(g)
disjoin(g)
coverage(g)

GRangesList的一个操作

先不学习，以后用到了再说

##################GenomicRanges的小练习##########################

读取玉米的注释文件

library(rtracklayer)
maize_gff<-import.gff(con = "E:\reference\Zea_mays.AGPv4.38.gff")
head(maize_gff)

提取出外显子区域

maize_gff_exon<-maize_gff[maize_gff$type=="exon",]
head(maize_gff_exon)

添加染色体长度信息

maize_gff_length<-maize_gff[maize_gff$type=="chromosome"]
maize_gff_length<-sort(maize_gff_length)
maize_gff_length
end(maize_gff_length)
seqlengths(maize_gff)<-end(maize_gff_length)

library(ggbio)
autoplot(seqinfo(maize_gff))

看一下外显子的长度和分布

width(maize_gff_exon)
hist(width(maize_gff_exon))
hist(log2(width(maize_gff_exon)))

################## ggbio画图 ###################
library(ggbio)
library(GenomicRanges)
install.packages("stringi")

p_ideo <- Ideogram(genome = "hg19")
p_ideo
p_ideo + xlim(GRanges("chr2", IRanges(1e8, 1e8+10000)))
Ideogram(genome = "hg19", xlabel = TRUE)

画详细的gene model图

画图时，数据可以来自OrganismDb，GRangesList，EnsDb，(TxDb)

1、使用organismdb的数据

library(Homo.sapiens)
class(Homo.sapiens)

data(genesymbol, package = "biovizBase")
wh <- genesymbol[c("BRCA1", "NBR1")]
wh <- range(wh, ignore.strand = TRUE)
p.txdb <- autoplot(Homo.sapiens, which = wh)
p.txdb
autoplot(Homo.sapiens, which = wh, label.color = "black", color = "brown",
fill = "brown")

使用gap.geom可以改变内含子的形状

autoplot(Homo.sapiens, which = wh, gap.geom = "chevron")
autoplot(Homo.sapiens, which = wh, stat = "reduce")
autoplot(Homo.sapiens, which = wh, columns = c("TXNAME", "GO"), names.expr = "TXNAME::GO")

2、使用txdb的数据

因为TxDb没有基因的symbol信息，但是可以根据别的数据的信息来画图

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
autoplot(txdb, which = wh)

3、使用ensdb的数据

source("http://bioconductor.org/biocLite.R")
biocLite("EnsDb.Hsapiens.v75")
library(EnsDb.Hsapiens.v75)
ensdb <- EnsDb.Hsapiens.v75

指定具体的基因

autoplot(ensdb, GenenameFilter("PHKG2"))
autoplot(ensdb, ~ symbol == "PHKG2", names.expr="gene_name")

指定区间

gr <- GRanges(seqnames = 16, IRanges(30768000, 30770000), strand = "*")
autoplot(ensdb, GRangesFilter(gr), names.expr = "gene_name")

使用基因id

autoplot(ensdb, GeneIdFilter(c("ENSG00000196118", "ENSG00000156873")))

4、使用GRangesList的数据

library(biovizBase)
gr.txdb <- crunch(txdb, which = wh)
colnames(values(gr.txdb))[4] <- "model"

grl <- split(gr.txdb, gr.txdb$tx_id)
names(grl) <- sample(LETTERS, size = length(grl), replace = TRUE)
grl

autoplot(grl, aes(type = model))
ggplot() + geom_alignment(grl, type = "model")

添加reference到track

library(BSgenome.Hsapiens.UCSC.hg19)
bg <- BSgenome.Hsapiens.UCSC.hg19
p.bg <- autoplot(bg, which = wh)

p.bg
p.bg + zoom(1/100)
p.bg + zoom(1/1000)
p.bg + zoom(1/2500)

autoplot(bg, which = resize(wh, width = width(wh)/2000), geom = "segment")

添加比对track

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- keepSeqlevels(wh, "chr17")
autoplot(fl.bam, which = wh)

显示gap

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- keepSeqlevels(wh, "chr17")
autoplot(fl.bam, which = resize(wh, width = width(wh)/10), geom = "gapped.pair")

显示错配

library(BSgenome.Hsapiens.UCSC.hg19)
bg <- BSgenome.Hsapiens.UCSC.hg19
p.mis <- autoplot(fl.bam, bsgenome = bg, which = wh, stat = "mismatch")
p.mis

显示覆盖情况

autoplot(fl.bam, method = "estimate")
autoplot(fl.bam, method = "estimate", which = paste0("chr", 17:18), aes(y = log(..coverage..)))

添加变异信息

可以使用variantanntoation来导入vcf文件，转换为vranges格式文件，

library(VariantAnnotation)
fl.vcf <- system.file("extdata", "17-1409-CEU-brca1.vcf.bgz", package="biovizBase")
vcf <- readVcf(fl.vcf, "hg19")

vr <- as(vcf[, 1:3], "VRanges")
vr <- renameSeqlevels(vr, value = c("17" = "chr17"))

gr17 <- GRanges("chr17", IRanges(41234400, 41234530))
p.vr <- autoplot(vr, which = wh)

p.vr
p.vr + xlim(gr17)
p.vr + xlim(gr17) + zoom()
autoplot(vr, which = wh, geom = "rect", arrow = FALSE)

gr17 <- GRanges("chr17", IRanges(41234415, 41234569))
tks <- tracks(p.ideo, mismatch = p.mis, dbSNP = p.vr, ref = p.bg, gene = p.txdb,
heights = c(2, 3, 3, 1, 4)) + xlim(gr17) + theme_tracks_sunset()
tks
tks + zoom()

p.txdb + zoom(1/8)
p.txdb + zoom(2)
p.txdb + nextView()
p.txdb + prevView()

画圆形图

画manhattan图

增加关于染色体图布局的数据

data(darned_hg19_subset500, package = "biovizBase")
dn <- darned_hg19_subset500
library(GenomicRanges)
seqlengths(dn)
head(dn)

data(ideoCyto, package = "biovizBase")
seqlengths(dn) <- seqlengths(ideoCyto$hg19)[names(seqlengths(dn))]
dn <- keepSeqlevels(dn, paste0("chr", c(1:22, "X")))
autoplot(dn, layout = "karyogram")

autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg))

autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg), alpha = 0.5) +
scale_color_discrete(na.value = "brown")

去除NA值

dn.nona <- dn[!is.na(dnlevels <- as.numeric(factor(dn.nona$exReg))
p.ylim <- autoplot(dn.nona, layout = "karyogram", aes(color = exReg, fill = exReg,
ymin = (levels - 1) * 10/3,
ymax = levels * 10 /3))
p.ylim
p.ylim + facet_wrap(~seqnames)

也可以使用layout_karyogram来添加layer

dn3 <- dn.nona[dn.nonaexReg == '5']
dnC <- dn.nona[dn.nonaexReg)]

autoplot(seqinfo(dn3), layout = "karyogram") +
layout_karyogram(data = dn3, geom = "rect", ylim = c(0, 10/3), color = "#7fc97f") +
layout_karyogram(data = dn5, geom = "rect", ylim = c(10/3, 10/32), color = "#beaed4") +
layout_karyogram(data = dnC, geom = "rect", ylim = c(10/32, 10), color = "#fdc086") +
layout_karyogram(data = dn.na, geom = "rect", ylim = c(10, 10/3*4), color = "brown")

dn$pvalue <- runif(length(dn)) * 10 #增加metadata
p <- autoplot(seqinfo(dn)) + layout_karyogram(dn, aes(x = start, y = pvalue),
geom = "point", color = "#fdc086")
p
p + theme_alignment()
p + theme_clear()
p + theme_null()

library(GenomicRanges)
set.seed(1)
N <- 100
gr <- GRanges(seqnames = sample(c("chr1", "chr2", "chr3"),
size = N, replace = TRUE),
IRanges(start = sample(1:300, size = N, replace = TRUE),
width = sample(70:75, size = N,replace = TRUE)),
strand = sample(c("+", "-"), size = N, replace = TRUE),
value = rnorm(N, 10, 3), score = rnorm(N, 100, 30),
sample = sample(c("Normal", "Tumor"),
size = N, replace = TRUE),
pair = sample(letters, size = N,
replace = TRUE))
seqlengths(gr) <- c(400, 1000, 500)
autoplot(gr)

autoplot(gr) + theme_genome()