annotate ChipSeqRatioAnalysis.R @ 12:6b7d31026d1c draft

Uploaded
author petr-novak
date Tue, 04 Feb 2020 08:38:22 -0500
parents 99569eccc583
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
1 #!/usr/bin/env Rscript
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
2 library(R2HTML, quietly=T)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
3 library(base64enc, quietly=T)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
4
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
5
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
6 htmlheader=
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
7 " <html xmlns:mml=\"http://www.w3.org/1998/Math/MathML\">
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
8 <head>
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
9 <title> ChIP-Seq Mapper Output </title>
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
10 <style>
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
11 html,body{font-family:Verdana,sans-serif;font-size:15px;line-height:1.5}
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
12
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
13 table {
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
14 border-collapse: collapse;
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
15 border: 1px solid black;
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
16 width: 1000pt
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
17 }
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
18 table, th, td {
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
19 border: 1px solid black;
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
20 }
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
21 </style>
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
22
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
23 </head>
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
24
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
25
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
26
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
27 "
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
28
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
29
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
30 #arguments
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
31 args <- commandArgs(trailingOnly = TRUE)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
32 input <- args[1]
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
33 HTMLfile <- args[2]
8
99569eccc583 Uploaded
petr-novak
parents: 6
diff changeset
34 thr = 5
99569eccc583 Uploaded
petr-novak
parents: 6
diff changeset
35 threshld <- thr/(thr+1)
99569eccc583 Uploaded
petr-novak
parents: 6
diff changeset
36
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
37 inputN=as.numeric(args[3])
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
38 chipN=as.numeric(args[4])
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
39 #dataframe preprocessing and table creation
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
40 df <- read.delim(input, comment.char="#")
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
41
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
42 df$"Ratio Chip/Input"=df$Chip_Hits/df$Input_Hits
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
43 df$"Normalized ratio Chip/Input"=(df$Chip_Hits/chipN)/(df$Input_Hits/inputN)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
44
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
45 df$"Ratio Chip/(Chip+Input)"=df$Chip_Hits/(df$Chip_Hits + df$Input_Hits)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
46 df$"Normalized ratio Chip/(Chip+Input)"=(df$Chip_Hits/chipN)/((df$Input_Hits/inputN)+(df$Chip_Hits/chipN))
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
47
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
48 outputTable = df[df$"Normalized ratio Chip/(Chip+Input)" > threshld,
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
49 ]
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
50 outputTable = outputTable[!is.na(outputTable$Cluster),
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
51 c('Cluster', 'Chip_Hits', 'Input_Hits',
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
52 'Normalized ratio Chip/Input','Normalized ratio Chip/(Chip+Input)')]
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
53 save.image("tmp.RData") #Plot creation
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
54 pngfile <- tempfile()
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
55 png(pngfile, width = 1000, height = 1200, pointsize=20)
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
56 par(mfrow=c(2,1))
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
57 lims=range(df$"Normalized ratio Chip/Input"[df$"Normalized ratio Chip/Input">0], finite = TRUE)
6
f224513123a1 Uploaded
petr-novak
parents: 5
diff changeset
58 suppressWarnings(plot(df$Cluster,df$"Normalized ratio Chip/Input", log="y", xlab="Cluster Nr.", ylab="Normalized ChiP/Input ratio", pch=20, ylim=lims))
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
59 abline(h=1,col='#00000080', lwd = 2)
8
99569eccc583 Uploaded
petr-novak
parents: 6
diff changeset
60 abline(h=thr,col='#FF000080', lwd = 2)
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
61
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
62
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
63 suppressWarnings(plot(df$Cluster,df$"Normalized ratio Chip/(Chip+Input)", xlab="Cluster Nr.", ylab="Normalized Chip/(Chip+Input)", pch=20))
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
64 abline(h=0.5,col='#00000080', lwd = 2)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
65 abline(h=threshld,col='#FF000080', lwd = 2)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
66
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
67
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
68 dev.off()
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
69 graph <- paste('<img src="data:image/png;base64 ,',
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
70 base64encode(pngfile),
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
71 '" alt="image" />'
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
72 )
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
73
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
74 #HMTL report creation + writing final output
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
75 directory=dirname(HTMLfile)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
76 filename=basename(HTMLfile)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
77 ## create HTML header
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
78 cat(htmlheader, file = filename)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
79
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
80
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
81 HTML(graph, file=filename)
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
82 if (nrow(outputTable)>0){
5
378565f5a875 Uploaded
petr-novak
parents: 3
diff changeset
83 HTML(outputTable, file=filename, classtable = "dataframe",
8
99569eccc583 Uploaded
petr-novak
parents: 6
diff changeset
84 row.names=FALSE, align='left', caption=paste("Clusters with Normalized ChIP/Input ratio >", thr), captionalign="top")
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
85 }
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
86 HTMLEndFile(filename)
6
f224513123a1 Uploaded
petr-novak
parents: 5
diff changeset
87 # file.rename(from=filename, to=HTMLfile)
f224513123a1 Uploaded
petr-novak
parents: 5
diff changeset
88 system(sprintf("cp -r ./%s %s", filename, HTMLfile))
3
e320ef2d105a Uploaded
petr-novak
parents:
diff changeset
89 write.table(df, file=input, sep="\t", row.names = FALSE)