|
12
|
1 ########################################################
|
|
|
2 #
|
|
|
3 # creation date : 09/08/16
|
|
|
4 # last modification : 16/08/16
|
|
|
5 # author : Dr Nicolas Beaume
|
|
|
6 # owner : IRRI
|
|
|
7 #
|
|
|
8 ########################################################
|
|
|
9 log <- file(paste(getwd(), "log_QC.txt", sep="/"), open = "wt")
|
|
|
10 sink(file = log, type="message")
|
|
|
11 #######################################################
|
|
|
12 dataStats.nbNA <- function(x) {
|
|
|
13 return(length(x[is.na(x)]))
|
|
|
14 }
|
|
|
15 dataStats.matrixSize <- function(data) {
|
|
|
16 nbMarker <- nrow(data)
|
|
|
17 nbSample <- ncol(data)
|
|
|
18 return(list(nbMarker=nbMarker, nbSample=nbSample))
|
|
|
19 }
|
|
|
20
|
|
|
21 parseFrq <- function(freqFile) {
|
|
|
22 data <- read.table(freqFile, h=T)
|
|
|
23 return(data.frame(chr=data$CHR, pos=data$SNP, maf=data$MAF))
|
|
|
24 }
|
|
|
25
|
|
|
26 parseHWE <- function(hweFile) {
|
|
|
27 data <- read.table(hweFile, h=T)
|
|
|
28 return(data.frame(chr=data$CHR, pos=data$SNP, pvalue=data$P))
|
|
|
29 }
|
|
|
30
|
|
|
31 ########################## main function ##################
|
|
|
32 createReport <- function(genoFile, freqFile, hweFile, out="report.txt") {
|
|
|
33 # get basic statistics (nb markers, nb samples, nb NA)
|
|
|
34 data <- read.table(genoFile, sep="\t", h=T)
|
|
|
35 dataDimension <- dataStats.matrixSize(data)
|
|
|
36 # get MAF info
|
|
|
37 freq <- parseFrq(freqFile)
|
|
|
38 # get Hardy-Weinberg equilibrium
|
|
|
39 hwe <- parseHWE(hweFile)
|
|
|
40 # merge info
|
|
|
41 info <- data.frame(freq, hwe=hwe[,3])
|
|
|
42 # write report
|
|
|
43 write(paste("number of marker :", dataDimension$nbMarker), file=out)
|
|
|
44 write(paste("number of sample :", dataDimension$nbSample), file=out, append = T)
|
|
|
45 write.table(info, file=out, append = T)
|
|
|
46 }
|
|
|
47 ############################ main ##########################
|
|
|
48 cmd <- commandArgs(T)
|
|
|
49 source(cmd[1])
|
|
|
50 createReport(genotype, plinkFreq, plinkHWE, out)
|
|
|
51 # # get basic statistics (nb markers, nb samples, nb NA)
|
|
|
52 # data <- read.table(genotype, sep="\t", h=T)
|
|
|
53 # dataDimension <- dataStats.matrixSize(data)
|
|
|
54 # # get MAF info
|
|
|
55 # freq <- parseFrq(freqFile)
|
|
|
56 # # get Hardy-Weinberg equilibrium
|
|
|
57 # hwe <- parseHWE(hweFile)
|
|
|
58 # # merge info
|
|
|
59 # info <- data.frame(freq, hwe=hwe[,3])
|
|
|
60 # # write report
|
|
|
61 # cat(paste("number of marker : ", dataDimension$nbMarker, "\n", sep=""))
|
|
|
62 # cat(paste("number of sample : ", dataDimension$nbSample, "\n", sep=""))
|
|
|
63 # print(info)
|
|
|
64 # cat("\n")
|
