Mercurial > repos > nicolas > oghma

comparison qualityControl.R @ 12:87f772c49a20 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author nicolas
date Fri, 21 Oct 2016 06:28:49 -0400
parents
children f818e787d0c0
comparison
equal deleted inserted replaced
11:502fbb316d2d 12:87f772c49a20
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")