沒空改了辽故，等有空了再來修改格式吧
############genomicranges的學(xué)習(xí)################

首先使用granges來創(chuàng)建一個對象

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))

查看對象中的內(nèi)容

seqnames(gr) #查看基因的名稱
ranges(gr) #查看基因的區(qū)域
strand(gr) #查看基因位于那條鏈
start(gr)
end(gr)
width(gr)
names(gr)
length(gr)

任何文件系統(tǒng)中的數(shù)據(jù)分為數(shù)據(jù)和元數(shù)據(jù)

數(shù)據(jù)是指普通文件中的實際數(shù)據(jù)物赶，

而元數(shù)據(jù)指用來描述一個文件的特征的系統(tǒng)數(shù)據(jù)魂毁，諸如訪問權(quán)限寿酌、文件擁有者以及文件數(shù)據(jù)塊的分布信息等等蒋得，

而在此處元數(shù)據(jù)就是指基因的附加信息

granges(gr) #查看基因基本信息
as.data.frame(granges(gr))#將基因的基本信息轉(zhuǎn)換為數(shù)據(jù)框

mcols(gr) #查看基因的元信息
as.data.frame(mcols(gr))#將基因的元數(shù)據(jù)轉(zhuǎn)換為數(shù)據(jù)框

如果要查看單獨的元數(shù)據(jù)乙墙，可以使用$

gr $score gr$ GC
mcols(gr) $score 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的元數(shù)據(jù)列

將某一個基因進行重復(fù)

singles<-split(gr,names(gr))#將gr根據(jù)names來分割
grMod<-gr[1:4]
grMod
grMod[2]<-singles[[1]]
grMod

grMod[2]<-grMod[1]

repeat, reverse, 重復(fù)和反向

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的方法，分為內(nèi)部咸产，外部，之間

內(nèi)部flank, resize, shift

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

如果移動過有負數(shù)脑溢，可以將其變?yōu)?

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

外部reduce，gaps赖欣，range屑彻，disjoin，coverage

reduce會對區(qū)間進行合并overlap畏鼓，得到一個簡化的區(qū)間

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

GRangesList的一個操作

先不學(xué)習(xí)酱酬，以后用到了再說

##################GenomicRanges的小練習(xí)##########################

讀取玉米的注釋文件

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

提取出外顯子區(qū)域

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圖

畫圖時，數(shù)據(jù)可以來自O(shè)rganismDb云矫，GRangesList，EnsDb汗菜，(TxDb)

1让禀、使用organismdb的數(shù)據(jù)

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可以改變內(nèi)含子的形狀

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的數(shù)據(jù)

因為TxDb沒有基因的symbol信息陨界，但是可以根據(jù)別的數(shù)據(jù)的信息來畫圖

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

3巡揍、使用ensdb的數(shù)據(jù)

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")

指定區(qū)間

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的數(shù)據(jù)

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來導(dǎo)入vcf文件菌瘪，轉(zhuǎn)換為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圖

增加關(guān)于染色體圖布局的數(shù)據(jù)

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(dn $exReg)] dn.nona$ levels <- 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.nona $exReg == '3'] dn5 <- dn.nona[dn.nona$ exReg == '5']
dnC <- dn.nona[dn.nona $exReg == 'C'] dn.na <- dn[is.na(dn$ exReg)]

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()