@@ -59,18 +59,36 @@ echo $CMD2
5959eval $CMD2
6060
6161
62- # annotate using R
62+ # # annotate using R
6363R --vanilla << RSCRIPT
6464options(stringsAsFactors=F)
6565#get data
6666dat.all <- read.delim("$InpFil ", colClasses="character")
6767
6868dat <- dat.all[,1:(grep("AlternateAlleles", colnames(dat.all)))]
69+ # Add a AAF column for Wen-I
70+ ADcols <- grep("AD$", colnames(dat))
71+ GTcols <- grep("GT$", colnames(dat))
72+ for(j in 1:length(ADcols)){
73+ dat <- cbind(dat, 0)
74+ ADcol <- ADcols[j]
75+ GTcol <- GTcols[j]
76+ nam <- gsub("AD", "AAF", colnames(dat)[ADcol])
77+ colnames(dat)[ncol(dat)] <- nam
78+ for(i in 1:nrow(dat)){
79+ AAs <- as.numeric(strsplit(dat[i,GTcol], "/")[[1]])
80+ AAs <- AAs[AAs!=0]+1
81+ ADs <- as.numeric(strsplit(dat[i,ADcol], ",")[[1]])
82+ tot <- sum(ADs)
83+ AAFs <- round(ADs[AAs]/tot,5)
84+ if(length(AAFs)==0) AAFs <- 0
85+ dat[i,nam] <- paste(AAFs, collapse=",")
86+ }
87+ }
6988#add leading space to genotypes to deal with Excel date format issue and read depths to stop auto correction to numbers
70- GTADcols <- grep("GT$|AD$", colnames(dat))
71- for(i in GTADcols){
72- dat[,i] <- paste(" ", dat[,i], sep="")
73- dat[,i] <- gsub(",", ", ", dat[,i])
89+ for(i in c(GTcols, ADcols)){
90+ dat[,i] <- paste(" ", dat[,i], sep="")
91+ dat[,i] <- gsub(",", ", ", dat[,i])
7492}
7593#get annotation table and extract relevant lines
7694annot <- read.delim("$AnnTab ")
0 commit comments