|
1
|
1 ########################################################################
|
|
|
2 # JAMMv1.0.7rev1 is a peak finder for joint analysis of NGS replicates.
|
|
|
3 # Copyright (C) 2014-2015 Mahmoud Ibrahim
|
|
|
4 #
|
|
|
5 # This program is free software: you can redistribute it and/or modify
|
|
|
6 # it under the terms of the GNU General Public License as published by
|
|
|
7 # the Free Software Foundation, either version 3 of the License, or
|
|
|
8 # (at your option) any later version.
|
|
|
9 #
|
|
|
10 # This program is distributed in the hope that it will be useful,
|
|
|
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
13 # GNU General Public License for more details.
|
|
|
14 #
|
|
|
15 # You should have received a copy of the GNU General Public License
|
|
|
16 # along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
17 #
|
|
|
18 # Contact: mahmoud.ibrahim@mdc-berlin.de
|
|
|
19 ########################################################################
|
|
0
|
20
|
|
|
21
|
|
|
22 # ==============================================
|
|
|
23 # Parsing Arguments (source: phantom SPP script)
|
|
|
24 # ==============================================
|
|
|
25 args = commandArgs(trailingOnly = TRUE) # Read Arguments from command line
|
|
|
26 #nargs = length(args) # number of arguments
|
|
|
27
|
|
|
28 #Set arguments to default values
|
|
|
29 infile = NA # input file
|
|
|
30 out = NA # output directory
|
|
|
31
|
|
|
32 #Parsing arguments and storing values
|
|
|
33 for (each.arg in args) {
|
|
|
34 #input bed file
|
|
|
35 if (grepl('^-infile=',each.arg)) {
|
|
|
36 arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
|
|
|
37 if (! is.na(arg.split[2]) ) {
|
|
|
38 infile <- arg.split[2]
|
|
|
39 } else {
|
|
|
40 stop('No input file')
|
|
|
41 }
|
|
|
42 }
|
|
|
43 if (grepl('^-out=',each.arg)) {
|
|
|
44 arg.split <- strsplit(each.arg,'=',fixed=TRUE)[[1]] # split on =
|
|
|
45 if (! is.na(arg.split[2]) ) {
|
|
|
46 out <- arg.split[2]
|
|
|
47 } else {
|
|
|
48 stop('No output directory')
|
|
|
49 }
|
|
|
50 }
|
|
|
51 }
|
|
|
52
|
|
|
53 #=======================> DONE!
|
|
|
54
|
|
|
55
|
|
|
56 # ===================
|
|
|
57 # Do the processing
|
|
|
58 # ===================
|
|
|
59
|
|
|
60 #Read in variables
|
|
|
61 if (file.exists(infile)) {
|
|
|
62
|
|
|
63 message("@", infile)
|
|
|
64 xcorrs = read.table(infile, header=FALSE)
|
|
|
65 shifts = as.numeric(xcorrs$V2)
|
|
|
66
|
|
|
67
|
|
|
68 #how many chromosomes
|
|
|
69 num = length(xcorrs$V1)
|
|
|
70
|
|
|
71 if (length(xcorrs$V2) > 0) {
|
|
|
72 if (length(xcorrs$V2) > 1) {
|
|
|
73 #shift stats
|
|
|
74 avg = mean(shifts)
|
|
|
75 std = sd(shifts)
|
|
|
76 mod = as.numeric(names(sort(-table(shifts)))[1])
|
|
|
77
|
|
|
78 #find out fragment length
|
|
|
79 fraglength = 100 #default
|
|
|
80 default = 1 #was default shift assumed
|
|
|
81 if (length(shifts) != 0) {
|
|
|
82 fraglength = ceiling(mean(shifts))
|
|
|
83 default = 0
|
|
|
84 }
|
|
|
85 } else {
|
|
|
86 avg = xcorrs$V2
|
|
|
87 std = xcorrs$V2
|
|
|
88 mod = xcorrs$V2
|
|
|
89 fraglength = xcorrs$V2
|
|
|
90 default = 0
|
|
|
91 }
|
|
|
92 } else {
|
|
|
93 fraglength = 100 #default
|
|
|
94 }
|
|
|
95 } else {
|
|
|
96 fraglength = 100
|
|
|
97 default = 1
|
|
|
98 }
|
|
|
99 #=======================> DONE!
|
|
|
100
|
|
|
101
|
|
|
102
|
|
|
103
|
|
|
104
|
|
|
105 # =====================
|
|
|
106 # Write result to File
|
|
|
107 # =====================
|
|
|
108 if (default == 0) {
|
|
|
109 writethis = paste0("Number of chromosomes used:", num, "\nAverage strand shift:", avg, "\nStandard Deviation:", std, "\nMode:", mod, "\nFragment Length:", fraglength)
|
|
|
110 write(writethis, file = paste0(out, "/xcorrsummary.txt"))
|
|
|
111
|
|
|
112 if (length(xcorrs$V2) > 1) {
|
|
|
113 i = pdf(file=paste0(out, "/xcorrsummary.pdf"))
|
|
|
114 i = plot(density(shifts), main="Strand Shifts Distribution over Chromosomes")
|
|
|
115 i = dev.off()
|
|
|
116 }
|
|
|
117 message("Fragment Length:", fraglength, "\n")
|
|
|
118 } else if (default == 1) {
|
|
|
119
|
|
|
120 writethis = paste0("Number of chromosomes used:", "0", "\nAverage strand shift:", "NA", "\nStandard Deviation:", "NA", "\nMode:", "NA", "\nFragment Length:", fraglength)
|
|
|
121 write(writethis, file = paste0(out, "/xcorrsummary.txt"))
|
|
|
122
|
|
|
123 message("Fragment length estimation failed, default fragment length assumed:", fraglength, "\n")
|
|
|
124 }
|
|
|
125 #=======================> DONE!
|
