0
$\begingroup$

I tried to find differences between delta G values in amino acids between 9 species. Here is a link to the table that I read in. (The header lists the species): https://www.dropbox.com/s/5l1kxytx94hhlsk/SPU_002551.txt?dl=0

In my program, I cut selected the highest and lowest 25% of the data to compare, though I am not sure which Anova to use. The repeated measures Anova assuming Sphericity seems to require 3 columns of data though I only have 2 and cannot add a 3rd position column because everything has been mixed up by selecting parts of the data. I cannot use the simple One-way Anova either because the variance at 25% of the data is not constant.

Is there any other type of Anova I could use to check if the delta G values between species are different while also having Greenhouse-Geisser and Huynh-Feldt Corrections?

Example of nonsensical output from Univariate Type III Repeated-Measures ANOVA Assuming Sphericity:

SS num Df   Error SS den Df        F    Pr(>F)    
(Intercept) 9.6914e+10      1 2656941741     88 3209.875 < 2.2e-16 ***
design      5.0190e+06      8   15829627    704   27.901 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Mauchly Tests for Sphericity

       Test statistic    p-value
design      0.0022698 4.6815e-87


Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity

       GG eps Pr(>F[GG])    
design 0.4192  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          HF eps   Pr(>F[HF])
design 0.4378048 6.190704e-18

These p-values are too small for the slight differences I am trying to detect

Code that produces the nonsensical result:

options(contrasts=c("contr.sum","contr.poly")) 
#adhere to the sum-to-zero convention for effect weights
#--------------------------------------------------------------------------------------------
#Getting input from anova.py
#--------------------------------------------------------------------------------------------
setwd("/Users/antonysagayaraj/Desktop/Anova")
info <- read.table("infoForRScript.txt", header=FALSE, sep = ',')

type <- toString(info[1,1])
cutLength <- as.numeric(info[1,2])
geneName <- toString(info[1,3])
setwd("/Users/antonysagayaraj/Desktop/Anova/CompiledText")
urchinframe <- read.table(paste(geneName,".txt",sep=""), header=TRUE, sep = ',')



#-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
#MAIN (A) Array
#-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

#-------------------------------------------------------------------------------------------
#Extracting columns as numeric + sorting 
#-------------------------------------------------------------------------------------------
sortAscending <- function(SpeciesArray) { #sorts the dG values from the bottom up
    return(sort(SpeciesArray, decreasing = FALSE, na.last = NA))
}

sortDescending <- function(SpeciesArray) { #sorts the dG values from the top down
    return(sort(SpeciesArray, decreasing = TRUE, na.last = NA))
}

cutAscending <- function(SpeciesArray,cutLength2) { #cuts the highest dG values
    temp <- sortAscending(SpeciesArray)
    return(temp[1:cutLength2])
}

cutDescending <- function(SpeciesArray,cutLength2) { #cuts the lowerst dG values
    temp <- sortDescending(SpeciesArray)
    return(temp[1:cutLength2])
}


Afragilis1 <- urchinframe[,c("A..fra")]
#print(Afragilis1)

Hpulcherrimus2 <- urchinframe[,c("H..pul")]
#print(Hpulcherrimus2)


Pdepressus3 <- urchinframe[,c("P..dep")]
#print(Pdepressus3)



Sdroebachiensis4 <- urchinframe[,c("S..dro")]
#print(Sdroebachiensis4)


Sfranciscanus5 <- urchinframe[,c("S..fran")]

#print(Sfranciscanus5)


Sintermedius6 <- urchinframe[,c("S..int")]
#print(Sintermedius6)


Snudus7 <- urchinframe[,c("S..nud")]
#print(Snudus7)



Spallidus8 <- urchinframe [,c("S..pal")]
#print(Spallidus8)


Spurpuratus9 <- urchinframe [,c("S..prp")]
#print(Spurpuratus9)


if (type == "percentage") { #makes sure that the percentage is converted into an amino acid number to be cut
    cutLength <- cutLength*length(Afragilis1)
}









#ANOVA FOR ASCENDING----------------------------------------------------------------------------------------Uses the top part of the data
dGlength <- length(cutAscending(Afragilis1,cutLength)) #gets the length of the cut data
ascendingaovurchin = data.frame(
    c(rep("Afra",dGlength),rep("Hpul",dGlength),rep("Pdep",dGlength),rep("Sdro",dGlength),rep("Sfra",dGlength),rep("Sint",dGlength),rep("Snud",dGlength),rep("Spal",dGlength),rep("Sprp",dGlength)))
    # makes a dataframe for use in the anova

colnames(ascendingaovurchin)= c("Species")
speciesList = list(cutAscending(Afragilis1,cutLength),cutAscending(Hpulcherrimus2,cutLength),cutAscending(Pdepressus3,cutLength),cutAscending(Sdroebachiensis4,cutLength),cutAscending(Sfranciscanus5,cutLength),cutAscending(Sintermedius6,cutLength),cutAscending(Snudus7,cutLength),cutAscending(Spallidus8,cutLength),cutAscending(Spurpuratus9,cutLength)) #list of sorted and cut dG values

for (x in 1:9) { #Goes down the species list and copies all the sorted and cut delta G values into the 2nd column of the data frame
  a <- x
  for (n in 1:(length(Sfranciscanus5))) {
    ascendingaovurchin[n + ((a-1)*dGlength),2] <- speciesList[[a]][n]
  }
}


matrix <- with(ascendingaovurchin, cbind(V2[Species=="Afra"], V2[Species=="Hpul"], V2[Species=="Pdep"], V2[Species=="Sdro"], V2[Species=="Sfra"], V2[Species=="Sint"], V2[Species=="Snud"], V2[Species=="Spal"], V2[Species=="Sprp"]))
model <- lm(matrix ~ 1)
design <- factor(c("Afra", "Hpul", "Pdep", "Sdro", "Sfra", "Sint", "Snud", "Spal", "Sprp"))
library(car)
aov <- Anova(model, idata=data.frame(design), idesign=~design, type="III")
summary(aov, multivariate=F)
|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ For each of your 9 columns you have several hundred rows. What do the rows represent? That's a lot more rows than there are amino acids. Restricting to the highest and lowest 25% is probably not a good way to proceed; unless there's a really good reason you should not go around throwing away data. In this case it would make intelligent interpretation of p-values and so forth almost impossible. $\endgroup$ – EdM Jul 17 '15 at 23:51
  • $\begingroup$ Yes. This text table represents one gene, which typically has several hundred amino acids. Each row is one position of a certain amino acid. We are trying to analyze where the differences between species are, and we need to restrict the data because only the peaks of the data are where the differences are. $\endgroup$ – anova Jul 18 '15 at 0:36
  • $\begingroup$ There presumably are differences in the amino acid sequences of this protein (not really a "gene") among your 9 species, or else the $\Delta G$ values would be the same at each position. Is there some reason why you are focused on the $\Delta G$ values instead of on the amino-acid sequences themselves? $\endgroup$ – EdM Jul 18 '15 at 2:08
  • $\begingroup$ It would help if you could be more specific about what the $\Delta G$ values represent. At first I though they might be $\Delta \Delta G$ values estimated theoretically for replacing each amino acid with, say, an alanine, but then there should be values of 0 at locations where there already are alanines. Without knowing what these $\Delta G$ values represent it will be hard to know if what you are trying to do makes sense and can actually be done. $\endgroup$ – EdM Jul 18 '15 at 9:48
1
$\begingroup$

"These $p$-values are too small for the slight differences I am trying to detect."

This statement is a little unusual because so long as you have enough data you can get $p$-values arbitrarily close to zero provided that there is some difference between the distributions. This is one criticism of the usefulness of $p$-values in general. As far as which ANOVA to use, I'm not entirely sure what you mean. It sounds like you have a fairly simple problem which can be addressed with any kind of two sample test. I would either run a $t$ test (which is of course essentially a special case of the one-way ANOVA) or Wilcoxon rank sum test. If the delta G values are multivariate in nature then you can use Hotelling's $T$-squared test, which generalizes the univariate $t$ test.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.