-
Notifications
You must be signed in to change notification settings - Fork 3
/
GAPIT.Genotype.R
executable file
·554 lines (420 loc) · 19 KB
/
GAPIT.Genotype.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
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
`GAPIT.Genotype` <-
function(G=NULL,GD=NULL,GM=NULL,KI=NULL,
kinship.algorithm=NULL,SNP.effect="Add",SNP.impute="Middle",PCA.total=0,seed=123, SNP.fraction =1,
file.path=NULL,file.from=NULL, file.to=NULL, file.total=NULL, file.fragment = 1000,SNP.test=TRUE,
file.G =NULL,file.Ext.G =NULL,
file.GD=NULL,file.Ext.GD=NULL,
file.GM=NULL,file.Ext.GM=NULL,
SNP.MAF=0.05,FDR.Rate = 0.05,SNP.FDR=1,
Timmer=NULL,Memory=NULL,
LD.chromosome=NULL,LD.location=NULL,LD.range=NULL, SNP.CV=NULL,
GP = NULL,GK = NULL,GTindex=NULL,
bin.size = 1000,inclosure.size = 100,
sangwich.top=NULL,sangwich.bottom=NULL,
file.output=TRUE,
Create.indicator = FALSE, Major.allele.zero = FALSE,Geno.View.output=TRUE){
#Object: To unify genotype and calculate kinship and PC if required:
# 1.For G data, convert it to GD and GI
# 2.For GD and GM data, nothing change
# 3.Samling GD and create KI and PC
# 4.Go through multiple files
# 5.In any case, GD must be returned (for QC)
#Output: GD, GI, GT, KI and PC
#Authors: Zhiwu Zhang
#Last update: August 11, 2011
##############################################################################################
#print("Genotyping: numericalization, sampling kinship, PCs and much more...")
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Genotype start")
Memory=GAPIT.Memory(Memory=Memory,Infor="Genotype start")
#Create logical variables
byData=!is.null(G) | !is.null(GD)
byFile=!is.null(file.G) | !is.null(file.GD)
hasGenotype=(byData | byFile )
needKinPC=(is.null(KI) | PCA.total>0 | kinship.algorithm=="Separation")
if(!is.null(KI) & !byData & !byFile & !SNP.test &kinship.algorithm!="SUPER") {
print("It return unexpected")
return (list(GD=NULL,GI=NULL,GT=NULL,hasGenotype=FALSE, genoFormat=NULL, KI=KI,PC=NULL,byFile=FALSE,fullGD=TRUE,Timmer=Timmer,Memory=Memory))
}
#Set indicator for full GD
fullGD=FALSE
if(byData) fullGD=TRUE
if(byFile & SNP.fraction==1 & needKinPC) fullGD=TRUE
#SET GT to NULL in case of no genotype
if(!byData & !byFile & is.null(GK) &kinship.algorithm!="SUPER") {
if(is.null(KI) & is.null(GP) & is.null(GK)) stop("GAPIT says: Kinship has to be provided or estimated from genotype!!!")
return (list(GD=NULL,GI=NULL,GT=NULL,hasGenotype=FALSE, genoFormat=NULL, KI=KI,PC=NULL,byFile=FALSE,fullGD=TRUE,Timmer=Timmer,Memory=Memory))
}
genoFormat="hapmap"
if(is.null(G)&is.null(file.G)) genoFormat="EMMA"
#Multiple genotype files
#In one of the 3 situations, calculate KI with the algorithm specified, otherwise skip cit by setting algorithm to "SUPER"
kinship.algorithm.save=kinship.algorithm
kinship.algorithm="SUPER"
#Normal
if(is.null(sangwich.top) & is.null(sangwich.bottom) ) kinship.algorithm=kinship.algorithm.save
#TOP or Bottom is MLM
pass.top=FALSE
if(!is.null(sangwich.top)) pass.top=!(sangwich.top=="FaST" | sangwich.top=="SUPER" | sangwich.top=="DC")
pass.bottom=FALSE
if(!is.null(sangwich.bottom)) pass.bottom=!(sangwich.bottom=="FaST" | sangwich.bottom=="SUPER" | sangwich.bottom=="DC")
if(pass.top | pass.bottom )kinship.algorithm=kinship.algorithm.save
#Compatibility of input
#agreement among file from, to and total
if(!is.null(file.from) &!is.null(file.to) &!is.null(file.total)){
if(file.total!=(file.to-file.from+1)) stop("GAPIT says: Conflict among file (from, to and total)")
}
if(!is.null(file.from) &!is.null(file.to)) {
if(file.to<file.from) stop("GAPIT says: file.from should smaller than file.to")
}
#file.from and file.to must be in pair
if(is.null(file.from) &!is.null(file.to) ) stop("GAPIT says: file.from and file.to must be in pair)")
if(!is.null(file.from) &is.null(file.to) ) stop("GAPIT says: file.from and file.to must be in pair)")
#assign file.total
if(!is.null(file.from) &!is.null(file.to) ) file.total=file.to-file.from+1
if(byFile& is.null(file.total)) stop("GAPIT says: file.from and file.to must be provided!)")
if(!is.null(GP) & !is.null(GK) ) stop("GAPIT Says: You can not provide GP and GK at same time")
if(!is.null(GP) & !is.null(KI) ) stop("GAPIT Says: You can not provide GP and KI at same time")
if(!is.null(GK) & !is.null(KI)) stop("GAPIT says: You can not specify GK and KI at same time!!!")
#GP does not allow TOP
if(!is.null(GP) & !is.null(sangwich.top) ) stop("GAPIT Says: You provided GP. You can not spycify sangwich.top")
#Top require a bottom
if(!is.null(sangwich.top) & is.null(sangwich.bottom) ) stop("GAPIT Says: Top require its Bottom")
#naked bottom require GP or GK
if(is.null(sangwich.top) & !is.null(sangwich.bottom) & (is.null(GP) & is.null(GK)) ) stop("GAPIT Says: Uncovered Bottom (without TOP) requires GP or GK")
#Pseudo top (GK or GP) requires a bottom
if(is.null(sangwich.top) & is.null(sangwich.bottom) & (!is.null(GP)|!is.null(GK ))) stop("GAPIT Says: You have provide GP or GK, you need to provide Bottom")
#if(!is.null(KI) &!is.null(kinship.algorithm)) stop("GAPIT says: You can not specify kinship.algorithm and provide kinship at same time!!!")
if(!needKinPC &SNP.fraction<1) stop("GAPIT says: You did not require calculate kinship or PCs. SNP.fraction should not be specified!!!")
if(!SNP.test & is.null(KI) & !byData & !byFile) stop("GAPIT says: For SNP.test optioin, please input either use KI or use genotype")
#if(is.null(file.path) & !byData & byFile) stop("GAPIT Ssays: A path for genotype data should be provided!")
if(is.null(file.total) & !byData & byFile) stop("GAPIT Ssays: Number of file should be provided: >=1")
if(!is.null(G) & !is.null(GD)) stop("GAPIT Ssays: Both hapmap and EMMA format exist, choose one only.")
if(!is.null(file.GD) & is.null(file.GM) & (!is.null(GP)|!is.null(GK)) ) stop("GAPIT Ssays: Genotype data and map files should be in pair")
if(is.null(file.GD) & !is.null(file.GM) & (!is.null(GP)|!is.null(GK)) ) stop("GAPIT Ssays: Genotype data and map files should be in pair")
if(!is.null(GD) & is.null(GM) & (is.null(GP)&is.null(GK)) &kinship.algorithm!="SUPER") stop("GAPIT Says: Genotype data and map files should be in pair")
if(is.null(GD) & !is.null(GM) & (is.null(GP)&is.null(GK)) &kinship.algorithm!="SUPER") stop("GAPIT Says: Genotype data and map files should be in pair")
#if(!byData & !byFile) stop("APIT Ssays: Either genotype data or files should be given!")
#if(byData&(!is.null(file.path))) stop ("APIT Ssays: You have provided geotype data. file.path should not be provided!")
#print("Pass compatibility of input")
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Genotype loaded")
Memory=GAPIT.Memory(Memory=Memory,Infor="Genotype loaded")
#Inital GLD
GLD=NULL
SNP.QTN=NULL #Intitial
GT=NULL
#Handler of read data in numeric format (EMMA)
#Rename GM as GI
if(!is.null(GM))GI=GM
rm(GM)
gc()
#Extract GD and GT from read data GD
if(!is.null(GD) )
{
GT=as.matrix(GD[,1]) #get taxa
GD=as.matrix(GD[,-1]) #remove taxa column
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="GT created from GD)")
Memory=GAPIT.Memory(Memory=Memory,Infor="GT created from GD")
}
#Hapmap format
if(!is.null(G))
{
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Before HapMap")
Memory=GAPIT.Memory(Memory=Memory,Infor="Before HapMap")
#Convert HapMap to numerical
print(paste("Converting genotype...",sep=""))
hm=GAPIT.HapMap(G,SNP.effect=SNP.effect,SNP.impute=SNP.impute, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero)
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="after HapMap")
Memory=GAPIT.Memory(Memory=Memory,Infor="after HapMap")
#Extracting SNP for LD plot
if(!is.null(LD.chromosome)){
#print("Extracting SNP for LD plot...")
chromosome=(G[,3]==LD.chromosome[1])
bp=as.numeric(as.vector(G[,4]))
deviation=abs(bp-as.numeric(as.vector(LD.location[1])) )
location=deviation< as.numeric(as.vector(LD.range[1]) )
index=chromosome&location
GLD=G[index,]
}else{
#print("No data in GLD")
GLD=NULL
}
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="HapMap")
Memory=GAPIT.Memory(Memory=Memory,Infor="HapMap")
print(paste("Converting genotype done.",sep=""))
#rm(G)
#gc()
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="G removed")
Memory=GAPIT.Memory(Memory=Memory,Infor="G removed")
GT=hm$GT
GD=hm$GD
GI=hm$GI
rm(hm)
gc()
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="hm removed")
Memory=GAPIT.Memory(Memory=Memory,Infor="hm removed")
}
#From files
if(!byData & byFile){
#print("Loading genotype from files...")
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="byFile")
Memory=GAPIT.Memory(Memory=Memory,Infor="byFile")
numFileUsed=file.to
if(!needKinPC)numFileUsed=file.from
#Initial GI as storage
GD=NULL
GT=NULL
GI=NULL
GLD=NULL
#multiple fragments or files
for (file in file.from:numFileUsed){
frag=1
numSNP=file.fragment
myFRG=NULL
#print(paste("numSNP before while is ",numSNP))
while(numSNP==file.fragment) { #this is problematic if the read end at the last line
print(paste("Reading file: ",file,"Fragment: ",frag))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Before Fragment")
Memory=GAPIT.Memory(Memory=Memory,Infor="Before Fragment")
myFRG=GAPIT.Fragment( file.path=file.path,file.from=file.from, file.to=file.to,file.total=file.total,file.G=file.G,file.Ext.G=file.Ext.G,
seed=seed,SNP.fraction=SNP.fraction,SNP.effect=SNP.effect,SNP.impute=SNP.impute,genoFormat=genoFormat,
file.GD=file.GD,file.Ext.GD=file.Ext.GD,file.GM=file.GM,file.Ext.GM=file.Ext.GM,
file.fragment=file.fragment,file=file,frag=frag,
LD.chromosome=LD.chromosome,LD.location=LD.location,LD.range=LD.range, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero)
#print(paste("numSNP after while is ",numSNP))
#print(paste("OK with file: ",file,"Fragment: ",frag))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="After Fragment")
Memory=GAPIT.Memory(Memory=Memory,Infor="After Fragment")
if(is.null(GT) & !is.null(myFRG$GT))GT= as.matrix(myFRG$GT)
if(is.null(GD)){
GD= myFRG$GD
}else{
if(!is.null(myFRG$GD)) {
GD=cbind(GD,myFRG$GD)
}
}
if(is.null(GI)){
GI= myFRG$GI
}else{
if(!is.null(myFRG$GI)) {
colnames(myFRG$GI)=c("SNP","Chromosome","Position")
GI=as.data.frame(rbind(as.matrix(GI),as.matrix(myFRG$GI)))
}
}
if(is.null(G)){
G= myFRG$G
}else{
if(!is.null(myFRG$G)) {
G=as.data.frame(rbind(as.matrix(G),as.matrix(myFRG$G[-1,])))
}
}
if(is.null(GLD)){
GLD= myFRG$GLD
}else{
if(!is.null(myFRG$GLD)) {
if(myFRG$heading){
GLD=as.data.frame(rbind(as.matrix(GLD),as.matrix(myFRG$GLD[-1,])))
}else{
GLD=as.data.frame(rbind(as.matrix(GLD),as.matrix(myFRG$GLD)))
}
}
}
#print("This fragment is joined")
if(file==file.from & frag==1)GT=as.matrix(myFRG$GT)
frag=frag+1
if(!is.null(myFRG$GI)) {
numSNP=myFRG$linesRead[1]
}else{
numSNP=0
}
if(!needKinPC)numSNP=0 #force to end the while loop
if(is.null(myFRG))numSNP=0 #force to end the while loop
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="END this Fragment")
Memory=GAPIT.Memory(Memory=Memory,Infor="END this Fragment")
} #end of repeat on fragment
# print("This file is OK")
} #end of file loop
print("All files loaded")
} #end of if(!byData&byFile)
#print("file loaded")
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Sampling genotype")
Memory=GAPIT.Memory(Memory=Memory,Infor="Sampling genotype")
#Plot thirt part kinship
if(!is.null(KI) &file.output) {
if(KI!=1) {
if(nrow(KI)<1000){
print("Plotting Kinship")
theKin=as.matrix(KI[,-1])
colnames(theKin)=KI[,1]
rownames(theKin)=KI[,1]
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="set kinship")
Memory=GAPIT.Memory(Memory=Memory,Infor="set kinship")
print("Creating heat map for kinship...")
pdf(paste("GAPIT.Kin.thirdPart.pdf",sep=""), width = 12, height = 12)
par(mar = c(25,25,25,25))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="prepare heatmap")
Memory=GAPIT.Memory(Memory=Memory,Infor="prepare heatmap")
heatmap.2(theKin, cexRow =.2, cexCol = 0.2, col=rev(heat.colors(256)), scale="none", symkey=FALSE, trace="none")
dev.off()
print("Kinship heat map PDF created!")
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="plot heatmap")
Memory=GAPIT.Memory(Memory=Memory,Infor="plot heatmap")
} #end of if(nrow(KI)<1000)
} #end of if(KI!=1)
} #end of if(!is.null(KI))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Before SUPER")
Memory=GAPIT.Memory(Memory=Memory,Infor="Before SUPER")
#SUPER
if(!is.null(GP) & kinship.algorithm=="SUPER" & !is.null(bin.size) & !is.null(inclosure.size)){
mySpecify=GAPIT.Specify(GI=GI,GP=GP,bin.size=bin.size,inclosure.size=inclosure.size)
SNP.QTN=mySpecify$index
if(!is.null(GD)){
#comment out to keep all taxa for GS, Zhiwu (Dec7, 2012)
#GK=GD[GTindex,SNP.QTN]
#SNPVar=apply(as.matrix(GK),2,var)
#GK=GK[,SNPVar>0]
#GK=cbind(as.data.frame(GT[GTindex]),as.data.frame(GK)) #add taxa
GK=GD[,SNP.QTN]
SNPVar=apply(as.matrix(GK),2,var)
GK=GK[,SNPVar>0]
GK=cbind(as.data.frame(GT),as.data.frame(GK)) #add taxa
#print("QTN extracted")
}
}
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Before PCA")
Memory=GAPIT.Memory(Memory=Memory,Infor="Before PCA")
#Create PC
PC=NULL
thePCA=NULL
if(PCA.total>0 | kinship.algorithm=="Separation"){
thePCA=GAPIT.PCA(X = GD, taxa = GT, PC.number = PCA.total,file.output=file.output,PCA.total=PCA.total)
PC=thePCA$PCs[,1:(1+PCA.total)]
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="PCA")
Memory=GAPIT.Memory(Memory=Memory,Infor="PCA")
print("PC created")
}
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Before creating kinship")
Memory=GAPIT.Memory(Memory=Memory,Infor="Before creating kinship")
#Create kinship from genotype if not provide
if(is.null(KI) & (!is.null(GD) |!is.null(GK)) & kinship.algorithm!="SUPER")
{
print("Calculating kinship...")
if(!is.null(GK)){
thisGD=GK[,-1]
myGT=as.matrix(GK[,1])
print("GK is used to create KI")
}else{
thisGD=GD
myGT=GT
#comment out to keep all taxa for GS, Zhiwu (Dec7, 2012)
#if(!is.null(GTindex)){
# thisGD=thisGD[GTindex,]
# myGT=myGT[GTindex]
#}
}
print(paste("Number of individuals and SNPs are ",nrow(thisGD)," and ",ncol(thisGD)))
theKin=NULL
if(kinship.algorithm=="EMMA"){
half.thisGD = as.matrix(.5*thisGD)
if(length(which(is.na(half.thisGD))) > 0){
print("Substituting missing values with heterozygote for kinship matrrix calculation....")
half.thisGD[which(is.na(half.thisGD))] = 1
}
theKin= emma.kinship(snps=t(as.matrix(.5*thisGD)), method="additive", use="all")
}
if(kinship.algorithm=="Loiselle")theKin= GAPIT.kinship.loiselle(snps=t(as.matrix(.5*thisGD)), method="additive", use="all")
if(kinship.algorithm=="VanRaden")theKin= GAPIT.kinship.VanRaden(snps=as.matrix(thisGD))
if(kinship.algorithm=="Zhang")theKin= GAPIT.kinship.ZHANG(snps=as.matrix(thisGD))
if(kinship.algorithm=="Separation")theKin= GAPIT.kinship.separation(PCs=thePCA$PCs,EV=thePCA$EV,nPCs=PCA.total)
if(!is.null(theKin)){
colnames(theKin)=myGT
rownames(theKin)=myGT
print("kinship calculated")
if(length(GT)<1000 &file.output){
#Create heat map for kinship
print("Creating heat map for kinship...")
pdf(paste("GAPIT.Kin.",kinship.algorithm,".pdf",sep=""), width = 12, height = 12)
par(mar = c(25,25,25,25))
heatmap.2(theKin, cexRow =.2, cexCol = 0.2, col=rev(heat.colors(256)), scale="none", symkey=FALSE, trace="none")
dev.off()
print("Kinship heat map created")
}
print("Adding IDs to kinship...")
#Write the kinship into a text file
KI=cbind(myGT,as.data.frame(theKin)) #This require big memory. Need a way to solve it.
print("Writing kinship to file...")
if(file.output) write.table(KI, paste("GAPIT.Kin.",kinship.algorithm,".csv",sep=""), quote = FALSE, sep = ",", row.names = FALSE,col.names = FALSE)
print("Kinship save as file")
rm(theKin)
gc()
}
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="Estimating kinship")
Memory=GAPIT.Memory(Memory=Memory,Infor="Estimating kinship")
print("Kinship created!")
} #end of if(is.null(KI)&!is.null(GD))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="after creating kinship")
Memory=GAPIT.Memory(Memory=Memory,Infor="after creating kinship")
#LD plot
#print("LD section")
if(!is.null(GLD) &file.output){
if(nrow(GLD)>500){
GLD=GLD[1,]
print("WARNING: The number of SNPs requested is beyond limitation. No LD plot created.")
}
if(nrow(GLD)>1)
{
print("Plot LD...")
hapmapgeno= data.frame(as.matrix(t(GLD[,-c(1:11)])))
hapmapgeno[hapmapgeno=="NN"]=NA
hapmapgeno[hapmapgeno=="XX"]=NA
hapmapgeno[hapmapgeno=="--"]=NA
hapmapgeno[hapmapgeno=="++"]=NA
hapmapgeno[hapmapgeno=="//"]=NA
LDdist=as.numeric(as.vector(GLD[,4]))
LDsnpName=GLD[,1]
colnames(hapmapgeno)=LDsnpName
#Prune SNM names
#LDsnpName=LDsnpName[GAPIT.Pruning(LDdist,DPP=7)]
LDsnpName=LDsnpName[c(1,length(LDsnpName))] #keep the first and last snp names only
#print(hapmapgeno)
print("Getting genotype object")
LDsnp=makeGenotypes(hapmapgeno,sep="",method=as.genotype) #This need to be converted to genotype object
print("Caling LDheatmap...")
pdf(paste("GAPIT.LD.chromosom",LD.chromosome,"(",round(max(0,LD.location-LD.range)/1000000),"_",round((LD.location+LD.range)/1000000),"Mb)",".pdf",sep=""), width = 12, height = 12)
#pdf(paste("GAPIT.LD.pdf",sep=""), width = 12, height = 12)
par(mar = c(25,25,25,25))
MyHeatmap <- try(LDheatmap(LDsnp, LDdist, LDmeasure="r", add.map=TRUE,
SNP.name=LDsnpName,color=rev(cm.colors(20)), name="myLDgrob", add.key=TRUE,geneMapLabelY=0.1) )
if(!inherits(MyHeatmap, "try-error")) {
#Modify the plot
grid.edit(gPath("myLDgrob", "Key", "title"), gp=gpar(cex=.5, col="blue")) #edit key title size and color
grid.edit(gPath("myLDgrob", "geneMap", "title"), gp=gpar(just=c("center","bottom"), cex=0.8, col="black")) #Edit gene map title
grid.edit(gPath("myLDgrob", "geneMap","SNPnames"), gp = gpar(cex=0.3,col="black")) #Edit SNP name
}else{
print("Warning: error in converting genotype. No LD plot!")
}
dev.off()
print("LD heatmap crated")
#grid.edit(gPath("myLDgrob", "heatMap","title"), gp=gpar(cex=1.0)) #Make title smaler
#grid.edit(gPath("myLDgrob", "geneMap", "title"), gp=gpar(just=c("right","bottom"), cex=0.5, col="blue")) #Edit gene map title
#grid.edit(gPath("myLDgrob", "Key", "labels"), gp=gpar(cex=.5, col="black")) #edit key lable size and color
}else{ # alternative of if(nrow(GLD)>1)
print("Warning: There are less than two SNPs on the region you sepcified. No LD plot!")
} #end of #if(nrow(GLD)>1)
}#end of if(!is.null(GLD))
Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="after LD plot")
Memory=GAPIT.Memory(Memory=Memory,Infor="after LD plot")
###output Marker density and decade of linkage disequilibrium over distance
if(!is.null(GI) & !is.null(GD) & file.output & Geno.View.output){
ViewGenotype<-GAPIT.Genotype.View(
myGI=GI,
myGD=GD,
#chr=1,
#w1_start=30,
#w1_end=230,
#mav1=10
)
}
#print("Genotype successfully acomplished")
return (list(G=G,GD=GD,GI=GI,GT=GT,hasGenotype=hasGenotype, genoFormat=genoFormat, KI=KI,PC=PC,byFile=byFile,fullGD=fullGD,Timmer=Timmer,Memory=Memory,SNP.QTN=SNP.QTN))
}
#=============================================================================================