I have a matrix (called results) that looks like this
id1 id2 id3 id4 id5 id6 id7 id8 id9
snp1 1 2 0 NA 1 1 1 2 1
snp2 2 2 2 2 0 2 NA NA 0
snp3 NA NA 1 NA 0 NA NA 2 2
So far, I have deleted rows and columns that were completely filled with NAs using
indexsnp=apply(results,1,
function(x) length(which(is.na(x)==T)))
indexsnp=which(indexsnp==length(results[1,]))
indexsample=apply(results,2,
function(x) length(which(is.na(x)==T)))
indexsample=which(indexsample==length(results[,1]))
#get rid of indexes
results=results[-indexsnp,]
results=results[,-indexsample]
I still have a lot of NAs in my dataset, so now I would like to see which snp have call rates below 95% (i.e. which rows consist of more than 5% NAs), and then delete those rows. I'm not sure how to do this. I have tried
snpsum.col <- col.summary(results)
library(snpStats)
call <- 0.95
use <- with(snpsum.col, (!is.na(Call.rate) & Call.rate >= call))
use[is.na(use)] <- FALSE
cat(ncol(results)-sum(use),"SNPs will be removed due to low call
rate.\n")
genotype <- genotype[,use]
snpsum.col <- snpsum.col[use,]
but I get the error
Error in col.summary(results) : not a SnpMatrix object
Is there another way I can do this?
If m is such a matrix, do
m <- m[is.na(m)%*%rep(1,ncol(m))<=ncol(m)*0.05,]