View Single Post
Old 03-10-2015, 05:35 AM   #3
youssefd
Junior Member
 
Location: Belgium

Join Date: Nov 2014
Posts: 4
Default

Hi,

I have analyzed the samples according to Devon's recommendation and everything went well without any error messages.

However, in the results, I do not understand why the interaction effect of Location and Status gives nearly the exact opposite log fold change.
Is this something normal or it indicates something wrong with the analysis?

Could it indicate a mislabeling of the samples? I am suspecting this already from the PCA plot of VST values.




Below is the relevant part of the code I used, starting from a SummarizedExperiment.

Code:
# Remove sample 4
se <- se[, colnames(se) != "SA4B4"]
se$SampleID <- factor(se$SampleID)
se$Patient <- factor(se$Patient)

# relevel: put NI and Ileum as base level
se$Status <- relevel(se$Status, "NI")
se$Location <- relevel(se$Location, "Ileum")

# create a DEseq dataset
dds <- DESeqDataSet(se = se, design = formula(~ Patient + Location + Status + Location:Status)) 
# run analysis
dds <- DESeq(dds)

# extract results 

# Test the hypothesis that here is no gene expression difference between Caecum and Ileum
results(dds, contrast=c("Location", "Caecum", "Ileum"))

# Test the hypothesis that here is no gene expression difference between I and NI 
results(dds, contrast=c("Status", "I", "NI"))

# Interaction effect between Caecum and Status
results(dds, contrast=list("LocationCaecum.StatusI", "LocationCaecum.StatusNI"))

# Interaction effect between Ileum and Status (The result of this contrast if the opposite of the contrast above)
results(dds, contrast=list("LocationIleum.StatusI", "LocationIleum.StatusNI"))

# Test the hypothesis that here is no gene expression difference between I and NI in Caecum
results(dds, contrast=list(c("StatusI", "LocationCaecum.StatusI"), c("StatusNI", "LocationCaecum.StatusNI")))

# Test the hypothesis that here is no gene expression difference between I and NI in Ileum
results(dds, contrast = list(c("StatusI", "LocationIleum.StatusI"), c("StatusNI", "LocationIleum.StatusNI")))
And below is the coefficients for the first 6 genes

Code:
head(coef(dds)[, 12:19])
DataFrame with 6 rows and 8 columns
                LocationIleum LocationCaecum    StatusNI     StatusI
                    <numeric>      <numeric>   <numeric>   <numeric>
ENSG00000000003  -0.137034259   0.1468754138 -0.06881193  0.07865309
ENSG00000000005  -0.483682079   0.4866318563 -0.08784116  0.09079094
ENSG00000000419   0.008710795   0.0003770367 -0.02530449  0.03439232
ENSG00000000457   0.121583429  -0.1132622163 -0.10059046  0.10891168
ENSG00000000460  -0.080019422   0.0865712307  0.11134169 -0.10478988
ENSG00000000938  -0.002435431   0.0101981947 -0.49175965  0.49952241
                LocationIleum.StatusNI LocationCaecum.StatusNI
                             <numeric>               <numeric>
ENSG00000000003            -0.11392008              0.11334093
ENSG00000000005            -0.42676013              0.42602082
ENSG00000000419            -0.09053654              0.09032357
ENSG00000000457             0.02507584             -0.02592245
ENSG00000000460            -0.14335772              0.14429482
ENSG00000000938             0.17662933             -0.18076818
                LocationIleum.StatusI LocationCaecum.StatusI
                            <numeric>              <numeric>
ENSG00000000003            0.11276674            -0.11210476
ENSG00000000005            0.42268926            -0.42192512
ENSG00000000419            0.09060986            -0.09032040
ENSG00000000457           -0.02405254             0.02496919
ENSG00000000460            0.14268424            -0.14356620
ENSG00000000938           -0.17664982             0.18085401
any help on this would be appreciated
Thanks in advance,
Youssef
youssefd is offline   Reply With Quote