I have RNA data generated using the smart-seq2 protocol and I am interested in running association analysis (some thing like logistic regression but with limma) between the gene expression values and cohort status. I have quantified the transcripts as TPM using RSEM. I am considering normalizing the data with TMM, estimating mean-variance relationship using voomWithQualityWeights and then using lmFit and eBayes to estimate the association.
I am not interested in seeing difference between cohort but to predict cohort status using gene expression values.
Following is my dummy code:
design <- model.matrix(~0+Cohort+sex+age, data=PhenoType) #Cohort:1/0 v <- voomWithQualityWeights(myNormalized_data, design=design, normalization="none", plot=TRUE) fit1 <- lmFit(v,design) fit2 <- eBayes(fit1)
All help and comments are highly appreciated.