File: qtltools_runFDR_atrans.R

package info (click to toggle)
qtltools 1.3.1%2Bdfsg-4
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, sid, trixie
  • size: 14,996 kB
  • sloc: cpp: 15,321; makefile: 243; sh: 63; ansic: 51
file content (26 lines) | stat: -rw-r--r-- 1,043 bytes parent folder | download | duplicates (2)
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
#!/usr/bin/env Rscript
#Load qvalue package
suppressMessages(library(qvalue))

#Read command line arguments
args <- commandArgs(trailingOnly = TRUE)
try(if(length(args) != 4) stop("Incorrect number of arguments, usage> Rscript runFDR_atrans.R adjusted.best.txt adjusted.hits.txt FDR output.txt"))

cat("\nProcessing QTLtools approximate trans output\n");
cat("  * File best  = [", args[1], "]\n");
cat("  * File hits  = [", args[2], "]\n");
cat("  * FDR        = [", args[3], "]\n");
cat("  * Output     = [", args[4], "]\n");


B = read.table(args[1], head=FALSE, stringsAsFactors=FALSE)
H = read.table(args[2], head=FALSE, stringsAsFactors=FALSE)
FDR = as.numeric(args[3])
B$qval = qvalue(B$V2)$qval
threshold = min(B$V2[which(B$qval > FDR)])
cat("  * Threshold of significance for adjusted P-values =" , threshold, "\n")

cat("\nFiltering hits and output results\n");
S = H[which(H$V8 <= threshold), ]
cat("  * " , nrow(S) , " are significante out of ", nrow(H), "\n")
write.table(S, args[4], quote=FALSE, row.names=FALSE, col.names=FALSE)