|
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 ########################################################
|
|
39
|
9
|
|
12
|
10 #######################################################
|
|
|
11 dataStats.nbNA <- function(x) {
|
|
|
12 return(length(x[is.na(x)]))
|
|
|
13 }
|
|
|
14 dataStats.matrixSize <- function(data) {
|
|
|
15 nbMarker <- nrow(data)
|
|
|
16 nbSample <- ncol(data)
|
|
|
17 return(list(nbMarker=nbMarker, nbSample=nbSample))
|
|
|
18 }
|
|
|
19
|
|
|
20 parseFrq <- function(freqFile) {
|
|
|
21 data <- read.table(freqFile, h=T)
|
|
|
22 return(data.frame(chr=data$CHR, pos=data$SNP, maf=data$MAF))
|
|
|
23 }
|
|
|
24
|
|
|
25 parseHWE <- function(hweFile) {
|
|
|
26 data <- read.table(hweFile, h=T)
|
|
|
27 return(data.frame(chr=data$CHR, pos=data$SNP, pvalue=data$P))
|
|
|
28 }
|
|
|
29
|
|
|
30 ########################## main function ##################
|
|
|
31 createReport <- function(genoFile, freqFile, hweFile, out="report.txt") {
|
|
|
32 # get basic statistics (nb markers, nb samples, nb NA)
|
|
|
33 data <- read.table(genoFile, sep="\t", h=T)
|
|
39
|
34 missingData <- dataStats.nbNA(data)/length(data)
|
|
12
|
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)
|
|
39
|
45 write(paste("ratio of missing data :", missingData), file=out, append = T)
|
|
|
46 suppressWarnings(suppressMessages(write.table(info, file=out, append = T)))
|
|
12
|
47 }
|
|
|
48 ############################ main ##########################
|
|
|
49 cmd <- commandArgs(T)
|
|
|
50 source(cmd[1])
|
|
|
51 createReport(genotype, plinkFreq, plinkHWE, out)
|
|
|
52 # # get basic statistics (nb markers, nb samples, nb NA)
|
|
|
53 # data <- read.table(genotype, sep="\t", h=T)
|
|
|
54 # dataDimension <- dataStats.matrixSize(data)
|
|
|
55 # # get MAF info
|
|
|
56 # freq <- parseFrq(freqFile)
|
|
|
57 # # get Hardy-Weinberg equilibrium
|
|
|
58 # hwe <- parseHWE(hweFile)
|
|
|
59 # # merge info
|
|
|
60 # info <- data.frame(freq, hwe=hwe[,3])
|
|
|
61 # # write report
|
|
|
62 # cat(paste("number of marker : ", dataDimension$nbMarker, "\n", sep=""))
|
|
|
63 # cat(paste("number of sample : ", dataDimension$nbSample, "\n", sep=""))
|
|
|
64 # print(info)
|
|
|
65 # cat("\n")
|
