Hello,
A colleague of mine brought to my attention that two GENESIS functions, fitNullMM() and assocTestMM() have become defunct in the latest version of GENESIS; and have been replaced by fitNullModel() and assocTestSingle(). I am now trying to work on modifying a crop GWAS pipeline that uses these defunct functions so that these new functions are incorporated. However, I am hitting a wall because of data formatting issues (for the genotype files) required for assocTestSingle().
I would report error message that I am getting, but I am not receiving any error messages. Instead, the object I am getting form assocTestSingle(iterator, nullmod, test = "Score") is empty (0 rows, and 9 columns).
I have deduced that the source of the error is likely from the formatting of the genotype file. That is, assocTestSingle() seems accept an object created by GenotypeBlockIterator() instead of an object created from GenotypeData() (and assocTestMM() ran fine for an object from GenotypeData). Therefore, my uneducated guess is that GenotypeBlockIterator() is not reading my data from GenotypeData().
Could I please trouble you to 30 minutes of your time over the next week or two for assistance with troubleshooting the formatting issues I am running into?
Thank you for your time.
Alex
Can you please share the code you are using, from creation of the GenotypeData object through the call to assocTestSingle?
Sure! Below is the code. Some quick notes: 0.) pheno.and.CV contains the phenotypes and covariates to include in the model. The first column is the indivdiual ID, the second column is the response variable, and then the remaining columns are covariates.
1.) hm.GD.without.snps.below.0.05.maf is an object where the rows are the SNPs and the columns are the SNP names. Each SNP is coded numerically (e.g., 0, 1, and 2 denote the homozygous, heterozygous, and homozygous genotypes, respectively)
2.) K is the kinship (= additive genetic relatedness) matrix. s
mydat.temp <- pheno.and.CV colnames(mydat.temp) <- c(“scanID”, “pheno”, “pc1”, “pc2”, “pc3”) mydat.temp$scanID <- as.character(mydat.temp$scanID) mydat.temp$pheno <- as.numeric(mydat.temp$pheno)
make ScanAnnotationDataFrame
scanAnnot <- ScanAnnotationDataFrame(mydat.temp) scanAnnot
Create the genotype file genotype <- t(hm.GD.without.snps.below.0.05.maf) snpID <- as.vector(as.integer(c(1:nrow(genotype)))) chromosome <- as.vectoras.integerthe.SNP.info[,2])) position <- as.vectoras.integerthe.SNP.info[,3])) scanID <- as.vector(as.integer(c(1:ncol(genotype))))
geno <- MatrixGenotypeReader(genotype = genotype, snpID = snpID, chromosome = chromosome, position = position, scanID = scanID) genoData <- GenotypeData(geno) iterator <- GenotypeBlockIterator(genoData, snpBlock = 1000000)
Take care of the kinship matrix mypcrel <- K[,-1] rownames(mypcrel) <- mydat.temp$scanID colnames(mypcrel) <- mydat.temp$scanID mypcrel <- as.matrix(mypcrel)
mypcrel.list <- list(“Kin” = mypcrel)
Fit the “null” model that does not include any tested markers
nullmod <- fitNullMM(scanData = scanAnnot, outcome = “pheno”, covars = c(“pc1”,”pc2”, “pc3”), covMatList = mypcrel.list,
family = binomial)
nullmod <- fitNullModel(scanAnnot, outcome = “pheno”, covars = c(“pc1”,”pc2”, “pc3”), cov.mat = mypcrel.list, family = binomial)
Problem: for some reason, the geno file will not accpet the line names. Solution: change the line names in the appropriate nullmod objects to the numeric ids
nullmod$scanID <- scanID rownames(nullmod$cholSigmaInv) <- scanID colnames(nullmod$cholSigmaInv) <- scanID
Run the GWAS scan
assoc <- assocTestMM(genoData = genoData, nullMMobj = nullmod, test = “Score”)
assoc <- assocTestSingle(iterator, nullmod, test = “Score”)
some Exploratory code
str(assoc) head(assoc)
the.results <- data.framethe.SNP.info$SNP, the.SNP.info$Chromosome, the.SNP.info$Position, assoc$Score.pval, assoc$Score.Stat)
P.S. - Sorry I am new to this bioconductor mailing list. Some of the code that I commented out (e.g., the now defunct assocTestMM() function) now appear in a big font.
If you edit your reply to use markdown format for your code blocks, it will be a lot easier to read: https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet#code-and-syntax-highlighting
Please elaborate on "Problem: for some reason, the geno file will not accept the line names." What error did you get exactly? Also, what version of GENESIS are you using?
Hi Stephanie,
Would you or somebody on your team please be available for a 30-minute maximum phone call for assistance some time over the next week? I think that the issues I am running into can be quickly and effectively resolved if I were to show my workspace, and the specific issues I am running into.
If this is possible, I would sincerely appreciate it.
Thank you again for your time and for developing such a great software package.
Best wishes,
Alex
I'm sorry, but I do not have time for phone consultations with package users. This support site exists for you to get help using Bioconductor packages, but you need to follow the guidelines. Please review the posting guidelines here: http://www.bioconductor.org/help/support/posting-guide/
Nevermind. I will figure this out myself.
Thank you for the time spent reading the email thread.