-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMakeDensityPlot.R
165 lines (135 loc) · 5.6 KB
/
MakeDensityPlot.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
library(karyoploteR)
library(RColorBrewer)
custom.genome <- toGRanges(data.frame(chr=c('C1','C2','C3','C4','C5','C6','C7','C8','C9'), start=c(1,1,1,1,1,1,1,1,1), end=c(43764888,52886895,64984695,53719093,46902585,39822476,48366697,41758685,54679868)))
gffGenes <- import.gff('Brassica_oleracea.v2.1.34.chr_with_PAV.gff3')
#head(gffGenes)
# These files were generated by RGAugury
gffNBS <- import.gff('Brassica.CViT.NBS.blue.txt')
gffRLK <- import.gff('Brassica.CViT.RLK.green.txt')
gffRLP <- import.gff('Brassica.CViT.RLP.orange.txt')
gffTMCC <- import.gff('Brassica.CViT.TMCC.black.txt')
#cols <- brewer.pal(6, 'Dark2')
gffGenes <- gffGenes[gffGenes$type == 'gene']
png('All_densitiesGenes.png', width=480*3, height=480*3)
kp <- plotKaryotype(genome = custom.genome, plot.type=4, ideogram.plotter = NULL, cex=1.4)
pp <- getDefaultPlotParams(plot.type=4)
pp$ideogramlateralmargin <- 0.009
kp <- plotKaryotype(genome = custom.genome, plot.type=4, cex=1.4, plot.params=pp)
kpAddBaseNumbers(kp, cex=1.4)
kpDataBackground(kp, data.panel = 1, r0=0, r1=0.2)
kp <- kpPlotDensity(kp, gffGenes, r0=0, r1=0.2, col=cols[1])#, window.size=10e6)
kpAddLabels(kp, labels='All\nGenes', pos=1, r0=0, r1=0.2, cex=1.4, label.margin=0.02)
kpDataBackground(kp, data.panel = 1, r0=0.23, r1=0.4)
kp <- kpPlotDensity(kp, gffNBS, r0=0.23, r1=0.4, col=cols[2])
kpAddLabels(kp, labels='NBS', pos=1, r0=0.23, r1=0.4, cex=1.4, label.margin=0.02)
kpDataBackground(kp, data.panel = 1, r0=0.43, r1=0.6)
kp <- kpPlotDensity(kp, gffRLK, r0=0.43, r1=0.6, col=cols[3])
kpAddLabels(kp, labels='RLK', pos=1, r0=0.43, r1=0.6, cex=1.4, label.margin=0.02)
kpDataBackground(kp, data.panel = 1, r0=0.63, r1=0.8)
kp <- kpPlotDensity(kp, gffRLP, r0=0.63, r1=0.8, col=cols[4])
kpAddLabels(kp, labels='RLP', pos=1, r0=0.63, r1=0.8, cex=1.4, label.margin=0.02)
kpDataBackground(kp, data.panel = 1, r0=0.73, r1=0.9)
kp <- kpPlotDensity(kp, gffTMCC, r0=0.73, r1=0.9, col=cols[5])
kpAddLabels(kp, labels='TMCC', pos=1, r0=0.73, r1=0.9, cex=1.4, label.margin=0.02)
sign.genes <- gffGenes[as.numeric(gffGenes$PAV)>1]
kp <- kpPoints(kp, data=sign.genes, y=as.numeric(sign.genes$PAV), ymin=0.0, ymax=100.0, r0=0.93, r1=1.1, col=cols[6])
kpDataBackground(kp, data.panel = 1, r0=0.83, r1=1)
kpAddLabels(kp, labels='Variable\ngenes', pos=1, r0=0.83, r1=1, cex=1.4, label.margin=0.02)
kp <- kpPlotDensity(kp, sign.genes, r0=0.83, r1=1, col=cols[6])
dev.off()
library(regioneR)
print('NBS')
pt <- permTest(A=sign.genes, B=gffNBS, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=numOverlaps, force.parallel=F, universe=gffGenes)
pt
png('NBS_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, A=sign.genes, B=gffNBS)
png('NBS_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
print(par('usr'))
dev.off()
print('RLK')
pt <- permTest(A=sign.genes, B=gffRLK, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=numOverlaps, force.parallel=F, universe=gffGenes)
pt
png('RLK_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, A=sign.genes, B=gffRLK)
png('RLK_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
print('RLP')
pt <- permTest(A=sign.genes, B=gffRLP, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=numOverlaps, force.parallel=F, universe=gffGenes)
pt
png('RLP_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, A=sign.genes, B=gffRLP)
png('RLP_lz.png')
plot(lz, pars=list(ylim=c(-1, 25)))
dev.off()
print('TMCC')
pt <- permTest(A=sign.genes, B=gffTMCC, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=numOverlaps, force.parallel=F, universe=gffGenes)
pt
png('TMCC_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, A=sign.genes, B=gffTMCC)
png('TMCC_lz.png')
plot(lz, pars=list(ylim=c(-1, 25)))
dev.off()
sign.genes <- sign.genes[sign.genes$type == 'gene']
te.genes <- import.gff('Order1_TO1000_Mito_Chloro_all_200bp_filtered_contaminants_marked.fasta.out_only_C.gff')
print(te.genes)
print('checking TE with NBS')
pt <- permTest(B=te.genes, A=gffNBS, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=meanDistance, universe=gffGenes)
pt
png('TE_NBS_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, B=te.genes, A=gffNBS)
png('NBS_TE_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
pt <- permTest(B=te.genes, A=gffRLK, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=meanDistance, universe=gffGenes, ntimes=500)
print('TE vs RLK')
pt
png('TE_RLK_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, B=te.genes, A=gffRLK)
png('RLK_TE_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
pt <- permTest(B=te.genes, A=gffRLP, genome=custom.genome, randomize.function=resampleRegions,evaluate.function=meanDistance, universe=gffGenes, ntimes=500)
print('TE vs RLP')
pt
png('TE_RLP_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, B=te.genes, A=gffRLP)
png('RLP_TE_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
pt <- permTest(B=te.genes, A=gffTMCC, genome=custom.genome, randomize.function=resampleRegions,evaluate.function=meanDistance, universe=gffGenes, ntimes=500)
print('TE vs TMCC')
pt
png('TE_TMCC_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, B=te.genes, A=gffTMCC)
png('TMCC_TE_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
print('checking TE with PAV')
pt <- permTest(B=te.genes, A=sign.genes, ntimes=500, genome=custom.genome, randomize.function=resampleRegions, evaluate.function=meanDistance, universe=gffGenes)
pt
png('TE_PAV_pt.png')
plot(pt)
dev.off()
lz <- localZScore(pt=pt, B=te.genes, A=gffNBS)
png('PAV_TE_lz.png')
plot(lz, pars=list(xlim=c(-1, 25)))
dev.off()
sessionInfo()