Question: I am having trouble recreating the linear regression results published in an Alzheimer's paper
0
gravatar for ravelarvargas
7 months ago by
ravelarvargas0 wrote:

I’m trying to recreate the results from this Alzheimer’s paper for further analysis (https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-016-0355-3)

Specifically, the supplementary table 4A reports genes significantly differentially expressed between brain donors that have Alzheimer’s vs. those that do not. The authors have corrected for covariates including sex, race, PMI, and pH for ~44,000 different probes (they did not remove half of the lowly expressed genes as is frequently suggested, as far as I’m aware).

I am having trouble recreating their results. Specifically, I am not getting any significant result when I look at the Inferior Frontal Gyrus differentially expressed genes when I look at low vs high Alzheimer’s clinical dementia rating (CDR). The problem is they are reporting 468 differentially expressed genes in the inferior frontal gyrus when they partition the data in the same way. I am hoping someone can help me pick up exactly what I am doing wrong.

I have two files, one is called patient and includes all of the patient covariates Patient covariates

The other file is the microarray expression data, called alzheimers Alzheimer's expression data

#The paper divided up patients into 'normal', 'low', and 'high' alzheimers (supplementary table 2 of paper) and I am dividing up the patients in the same way
CDR_low=patient%>%dplyr::filter(CDR%in%c(0.5, 1, 2))%>%mutate(alzheimers='low') 
CDR_high=patient%>%dplyr::filter(CDR%in%c(3, 4, 5))%>%mutate(alzheimers='high')

#bind patient tables together
CDR_low_high=rbind(CDR_low, CDR_high)

#filter out alzheimers table results that don't match the patient data
CDR_low_high_alzheimers=c()
CDR_low_high_alzheimers_colnames=c()
for (i in 1:ncol(alzheimers)){
  if (colnames(alzheimers)[i]%in%CDR_low_high$individualIdentifier){
    CDR_low_high_alzheimers=cbind(CDR_low_high_alzheimers, alzheimers[,i])
    CDR_low_high_alzheimers_colnames=rbind(CDR_low_high_alzheimers_colnames, colnames(alzheimers)[i])

  }
}

colnames(CDR_low_high_alzheimers)=CDR_low_high_alzheimers_colnames
rownames(CDR_low_high_alzheimers)=alzheimers$ID

#filter out patient table results that don't match the alzheimers data
CDR_low_high_present=c()
for (i in 1:nrow(CDR_low_high)){
  if (CDR_low_high$individualIdentifier[i]%in%colnames(CDR_low_high_alzheimers)){
    CDR_low_high_present=rbind(CDR_low_high_present, CDR_low_high[i,])
  }
}
CDR_low_high=CDR_low_high_present

rownames(CDR_low_high)=c()

rownames(CDR_low_high)=CDR_low_high$individualIdentifier
CDR_low_high=CDR_low_high%>%dplyr::select(PMI, pH, Sex, Race, Age, alzheimers)

# take log2 since the data is not normally distributed
CDR_low_high_alzheimers=log2(CDR_low_high_alzheimers)
rownames(CDR_low_high_alzheimers)=alzheimers$ID

#Design matrix
CDR_low_high_design <- model.matrix(~0+CDR_low_high$alzheimers+CDR_low_high$pH+CDR_low_high$Sex+CDR_low_high$Race+CDR_low_high$PMI)

colnames(CDR_low_high_design)=c('CDR_alzheimers_high', 'CDR_alzheimers_low', 'pH', 'sex_male', 'race_hispanic', 'race_white', 'PMI')

CDR_low_high_expression_fit <- lmFit(CDR_low_high_alzheimers, CDR_low_high_design)

CDR_low_high_expression_contr.matrix <- makeContrasts(
  CDR_alzheimers_low_CDR_alzheimers_normal = CDR_alzheimers_high-CDR_alzheimers_low, 
  levels = CDR_low_high_design)

CDR_low_high_expression_fit2 <- contrasts.fit(CDR_low_high_expression_fit, CDR_low_high_expression_contr.matrix)
CDR_low_high_expression_fit2 <- eBayes(CDR_low_high_expression_fit2)
topTable(CDR_low_high_expression_fit2, number=nrow(CDR_low_high_expression_fit2))

This is what my results look like: my differential gene expression results

this is what their results look like: their differential gene expression results

as you can see, they have a lot more significant results whereas I have none. Furthermore, even their log2FC values are different and I don't know what I am doing wrong or how I can replicate their results.

ADD COMMENTlink modified 7 months ago • written 7 months ago by ravelarvargas0

1) Did you contact the authors? They should provide code that you could study/use/improve.

2) You do not describe how you retrieve the data to run your script. Your script is long and hard to digest. I used GEOquery to obtain some of the data from GEO, and was able to do a sanity check of the assertion on TXNIP, which I think is the gene asserted to exhibit largest difference in the subset you are looking at.

ADD REPLYlink written 7 months ago by Vincent J. Carey, Jr.6.3k

1) I contacted four of the authors but unfortunately they all directed me to one of the four authors I contacted who hasn't gotten back to me over a few weeks now. I don't know if he is on holiday or busy but I am at a standstill until I can figure out how they got their results.

2) I got the data from the data portal on synapse, which you need an account to access. Specifically I am using this file (https://www.synapse.org/#!Synapse:syn16779991). the patient file with all the covariates is here (https://www.synapse.org/#!Synapse:syn3205399). The data itself has already been processed from the CEL files.

ADD REPLYlink written 7 months ago by ravelarvargas0
1

I will not be looking at synapse, so cannot help with that. You might want to work from the GEO deposit, checking that it is consistent with the synapse content. If the datasets are consistent, you would be able to post a fully reproducible example based on the GEO data and then you might get more assistance here. This is a bit of how I started out.

library(GEOquery)
gg = getGEO("GSE84422")
table(gg[[2]][["brain region:ch1"]])

            Anterior Cingulate                Caudate Nucleus 
                            59                             52 
Dorsolateral Prefrontal Cortex                   Frontal Pole 
                            57                             63 
                   Hippocampus         Inferior Frontal Gyrus 
                            55                             53 
       Inferior Temporal Gyrus          Middle Temporal Gyrus 
                            58                             58 
       Occipital Visual Cortex          Parahippocampal Gyrus 
                            53                             60 
    Posterior Cingulate Cortex               Precentral Gyrus 
                            58                             49 
             Prefrontal Cortex                        Putamen 
                            56                             52 
      Superior Parietal Lobule        Superior Temporal Gyrus 
                            50                             60 
                 Temporal Pole 
                            58 

ADD REPLYlink written 7 months ago by Vincent J. Carey, Jr.6.3k

I've added a fully reproducible version here https://support.bioconductor.org/p/118615/.

Still no luck with the original question unfortunately

ADD REPLYlink written 7 months ago by ravelarvargas0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 176 users visited in the last hour