Bioscripts/Plot_manhat_CircularDouble_CMplot.R at main · PotapenkoEugene/Bioscripts · GitHub

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
if(!require('dplyr')) install.packages('dplyr'); library('dplyr')
if(!require('data.table')) install.packages('data.table'); library('data.table')
source("https://raw.githubusercontent.com/YinLiLin/CMplot/master/R/CMplot.r")

args = commandArgs(trailingOnly=TRUE)
################################
GWAS1 = args[1] # LFMM gwasUniq custom format
GWAS2 = args[2] # EMMAX gwasUniq custom format
TRAIT = args[3]
SORTED = args[4] # T or F, if FALSE it will be sorted by script
if(length(args) == 5){
	CUSTOM_SNP = fread(args[5], header = F)$V1 # should be list without header
}else(CUSTOM_SNP = NULL)
################################


# Prepare DFs

message('INFO: PREPARING GWAS1 DATA FRAME')
lfmm.res = fread(GWAS1, sep = '\t', header = T) %>%
	setNames(c('SNP', 'Chromosome', 'Position', 'LFMM'))

message('INFO: PREPARING GWAS2 DATA FRAME')
emmax.res = fread(GWAS2, sep = '\t', header = T)

lfmm.res <-
	lfmm.res %>%
		left_join(emmax.res %>% dplyr::select(SNP, pvalue), by = 'SNP') %>%
		dplyr::rename(EMMAX = pvalue)

if(SORTED == 'F'){
	lfmm.res <-
		lfmm.res %>%
		dplyr::arrange(Chromosome)
}

if(is.null(CUSTOM_SNP)){
	pval.threshold = 0.05 / nrow(lfmm.res)

	# Plot
	message('INFO: PLOTTING WITH PVALUE THRESHOLD')
	CMplot(lfmm.res  %>% as.data.frame, # return some error to data.table
	       type="p",
	       #plot.type="c", #circular
	       #r=0.4,

	       plot.type='m',
	       multracks= T,

	       col=c("dodgerblue4", "deepskyblue"),
	       cex = 0.2,
	       signal.cex = 0.4,

	       main = TRAIT,

	       threshold= pval.threshold,
	       cir.chr.h=1.5,
	       amplify=T,
	       threshold.lty=1,
	       threshold.col="black",
	       signal.line=1,
	       signal.col="red",
	       chr.den.col=c("darkgreen","yellow","red"),
	       bin.size=1e6,
	       outward=FALSE,
	       file.name = paste0('CMplot_CircularDouble_EMMAXandLFMM_', TRAIT),
	       file="png",
	       dpi=600,
	       file.output=TRUE,
	       verbose=TRUE,

	       chr.labels.angle=60,
	       cex = 0.25,
	       axis.cex = 0.85,

	       width=10,height=5)
 }else{

	 # Plot
        message('INFO: PLOTTING WITH CUSTOM SNP HIGHLIGHTER')
        CMplot(lfmm.res  %>% as.data.frame, # return some error to data.table
               type="p",
               plot.type="c",
               r=0.4,
               col=c("dodgerblue4", "deepskyblue"),
               cex = 0.2,
               signal.cex = 0.4,

               #threshold= pval.threshold,
		highlight = CUSTOM_SNP,
	       	highlight.cex = 0.5,
		highlight.col = 'red',


               cir.chr.h=1.5,
               amplify=T,
               threshold.lty=1,
               threshold.col="black",
               signal.line=1,
               signal.col="red",
               chr.den.col=c("darkgreen","yellow","red"),
               bin.size=1e6,
               outward=FALSE,
               file.name = paste0('CMplot_CircularDouble_EMMAXandLFMM_customSNP_', TRAIT, '.png'),
	       file="png",
               dpi=1200,
               file.output=TRUE,
               verbose=TRUE,
               width=10,height=5)
}