Hi , im trying to preform Cochran–Armitage test on SNP data of 21 patients to find assosiation between disease sevirity( HIGH,MILD) and snp. for every snp im preforming the test and my gouls is to get for every sno the pvalue . for this im using the CATT function in R (https://search.r-project.org/CRAN/refmans/CATT/html/CATT.html) or this one :
catt <-
function(y, x, score = c(0, 1, 2)) {
miss <- unique(c(which(is.na(y)), which(is.na(x))))
n.miss <- length(miss)
if(n.miss > 0) {
y <- y[-miss]
x <- x[-miss]
}
if(!all((y == 0) | (y == 1)))
stop("y should be only 0 or 1.")
if(!all((x == 0) | (x == 1) |(x == 2)))
stop("x should be only 0, 1 or 2.")
ca <- x [y == 1]
co <- x [y == 0]
htca <- table(ca)
htco <- table(co)
A <- matrix(0, 2, 3)
colnames(A) <- c(0, 1, 2)
rownames(A) <- c(0, 1)
A[1, names(htca)] <- htca
A[2, names(htco)] <- htco
ptt <- prop.trend.test(A[1, ], colSums(A), score = score)
res=list("2x3-table" = A,
#chisq = as.numeric(ptt$statistic),
#df = as.numeric(ptt$parameter),
p.value=as.numeric(ptt$p.value))
#n.miss = n.miss)
return(res)
}
both result are the same . in my data i have an encoding column for the SNP (hetro=1,homo=2) when my data looks like this for example :
is_severe encoding
0 2
1 1
it gives me the pvalue:
The Cochran-Armitage Trend Test
data: The type of data is variable!
Z = 1.414, p-value = 0.1573
$`2x3-table`
0 1 2
0 0 1 0
1 0 0 1
but when my data looks like this:
is_severe encoding
0 2
0 1
print(catt(sn$ï..is_severe,sn$encoding))
Error in lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, :
NA/NaN/Inf in 'x'
> print(CATT(sn$ï..is_severe,sn$encoding))
Error in `[.default`(tbl, 2, ) : subscript out of bounds
also when data looks like this
is_severe encoding
0 1
0 1
it returns the same result
print(catt(sn$ï..is_severe,sn$encoding))
Error in lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, :
NA/NaN/Inf in 'x'
> print(CATT(sn$ï..is_severe,sn$encoding))
Error in `[.default`(tbl, 2, ) : subscript out of bounds
>
im new to genatics and i dont know how to fix this or what should i do with my data in order for the test to work on all SNPs and not part of them . thank you:)