Mercurial > repos > nicolas > oghma

comparison qualityControl.R @ 39:f818e787d0c0 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author nicolas
date Tue, 25 Oct 2016 14:41:59 -0400
parents 87f772c49a20
children
comparison
equal deleted inserted replaced
38:8112bc642858 39:f818e787d0c0
4 # last modification : 16/08/16 4 # last modification : 16/08/16
5 # author : Dr Nicolas Beaume 5 # author : Dr Nicolas Beaume
6 # owner : IRRI 6 # owner : IRRI
7 # 7 #
8 ######################################################## 8 ########################################################
9 log <- file(paste(getwd(), "log_QC.txt", sep="/"), open = "wt") 9
10 sink(file = log, type="message")
11 ####################################################### 10 #######################################################
12 dataStats.nbNA <- function(x) { 11 dataStats.nbNA <- function(x) {
13 return(length(x[is.na(x)])) 12 return(length(x[is.na(x)]))
14 } 13 }
15 dataStats.matrixSize <- function(data) { 14 dataStats.matrixSize <- function(data) {
30 29
31 ########################## main function ################## 30 ########################## main function ##################
32 createReport <- function(genoFile, freqFile, hweFile, out="report.txt") { 31 createReport <- function(genoFile, freqFile, hweFile, out="report.txt") {
33 # get basic statistics (nb markers, nb samples, nb NA) 32 # get basic statistics (nb markers, nb samples, nb NA)
34 data <- read.table(genoFile, sep="\t", h=T) 33 data <- read.table(genoFile, sep="\t", h=T)
34 missingData <- dataStats.nbNA(data)/length(data)
35 dataDimension <- dataStats.matrixSize(data) 35 dataDimension <- dataStats.matrixSize(data)
36 # get MAF info 36 # get MAF info
37 freq <- parseFrq(freqFile) 37 freq <- parseFrq(freqFile)
38 # get Hardy-Weinberg equilibrium 38 # get Hardy-Weinberg equilibrium
39 hwe <- parseHWE(hweFile) 39 hwe <- parseHWE(hweFile)
40 # merge info 40 # merge info
41 info <- data.frame(freq, hwe=hwe[,3]) 41 info <- data.frame(freq, hwe=hwe[,3])
42 # write report 42 # write report
43 write(paste("number of marker :", dataDimension$nbMarker), file=out) 43 write(paste("number of marker :", dataDimension$nbMarker), file=out)
44 write(paste("number of sample :", dataDimension$nbSample), file=out, append = T) 44 write(paste("number of sample :", dataDimension$nbSample), file=out, append = T)
45 write.table(info, file=out, append = T) 45 write(paste("ratio of missing data :", missingData), file=out, append = T)
46 suppressWarnings(suppressMessages(write.table(info, file=out, append = T)))
46 } 47 }
47 ############################ main ########################## 48 ############################ main ##########################
48 cmd <- commandArgs(T) 49 cmd <- commandArgs(T)
49 source(cmd[1]) 50 source(cmd[1])
50 createReport(genotype, plinkFreq, plinkHWE, out) 51 createReport(genotype, plinkFreq, plinkHWE, out)