Pair wise comparision using limma

Question

[iti.gupta]

Hello all, I know there have been a lot of post about pair-wise comparison using limma but still I am struggling for the same. Please help!!! I am facing issue in defining the design matrix and the downstream steps.

My experimental design looks like the following:

Patient    Tissue (same organ)

VD16         L

VD16         NL

VD24         L

VD24         NL

The codes so far:

tissue <-factor(targets$Tissue)

patient<-factor(targets$Sample_Name)

test_design<-factor(paste(patient, tissue, sep="."))

design <- model.matrix(~0+test_design)

[Benn]

I think your design should be:

patient <- factor(targets$Patient)

design <- model.matrix(~0+tissue+patient)

[iti.gupta]

My colnames for the design matrix looks like:

colnames(design)
 [1] "tissueL"     "tissueNL"    "patientVD18" "patientVD19" "patientVD20" "patientVD22"
 [7] "patientVD23" "patientVD24" "patientVD32" "patientVD33" "patientVD37" "patientVD38"
[13] "patientVD39" "patientVD40" "patientVD44" "patientVD49"

The downstream code will be like this:

cm<-makeContrasts(tissueL-tissueNL,levels = design)

fit<-lmFit(data_expressed,design)

fit1<-contrasts.fit(fit, contrasts = cm)

fit2<-eBayes(fit)

Is it correct? Also, how can I get the list of differentially expressed genes (consensus across all samples)?

[russhh]

Use topTable or decideTests to get the differnetially expressed genes. Have you modified the encoding of your tissue names? According to your original code there should be a tissueNormal and a tissueTumor column in the design matrix

[iti.gupta]

Sorry for the confusion I have corrected the tissue encoding in the above posts. Also, I have used the topTable But the resulting matrix is somewhat unusual.

colnames(tab)
 [1] "TargetID"    "tissueL"     "tissueNL"    "patientVD18" "patientVD19" "patientVD20"
 [7] "patientVD22" "patientVD23" "patientVD24" "patientVD32" "patientVD33" "patientVD37"
[13] "patientVD38" "patientVD39" "patientVD40" "patientVD44" "patientVD49" "AveExpr"    
[19] "F"           "P.Value"     "adj.P.Val"

I want a single column for logFC, P-value etc. for all samples combined(Somewhat like this: "logFC" "AveExpr" "t" "P.Value" "adj.P.Val" "B" )

[Benn]

Please follow the user's guide first instead of asking here if each step you take is correct. No offense. When you encounter problems, come back here and explain exactly what error you encounter and give use your exact code. Good luck.

[iti.gupta]

Hey ben, I have checked the user guide. I was able to compare the simple case versus control. However, I was facing the issue with paired comparisons hence i posted here. The complete code is posted below

      > targets<-readTargets("Targets.txt")
    > head(targets)
      Sample_ID Sample_Name Gender Tissue Age Age_Onset Merge_data
    1     VD13N        VD13      M     NL  21        18 VD13N:M:21
    2     VD13V        VD13      M      L  21        18 VD13V:M:21
    3     VD18N        VD18      F     NL  NA        NA VD18N:F:NA
    4     VD18V        VD18      F      L  NA        NA VD18V:F:NA
    5     VD19N        VD19      M     NL  21        NA VD19N:M:21
    6     VD19V        VD19      M      L  21        NA VD19V:M:21
    > data<-read.ilmn(files = "Probe_Profile.txt",ctrlfiles = "Control_Probe_Profile.txt", probeid = "ProbeID" , annotation = "TargetID", other.columns = c("Detection_Pval", "Avg_NBEADS"))
    Reading file Probe_Profile.txt ... ...
    Reading file Control_Probe_Profile.txt ... ...
    > dim(data)
    [1] 49587    30
    > data$E[1:5,1:5]
           VD13N    VD13V    VD18N    VD18V    VD19N
6450255 59.08015 61.36931 55.01596 64.49313 57.28073
2570615 78.69041 85.11665 69.07433 79.48322 80.94968
6370619 60.95150 59.82449 57.76508 59.30838 63.49848
2600039 68.59216 60.82525 60.31023 66.60255 63.30713
2650615 66.20662 69.16654 63.56881 63.85797 66.63638
    > table(data$genes$Status)

                BIOTIN            CY3_HYB       HOUSEKEEPING           LABELING 
                     2                  6                  7                  2 
    LOW_STRINGENCY_HYB           NEGATIVE            regular 
                     8                759              48803 
    > pe <- propexpr(data)
    > pe
        VD13N     VD13V     VD18N     VD18V     VD19N     VD19V     VD20N     VD20V     VD23N 
    0.5228984 0.5649895 0.5429903 0.5442841 0.5053918 0.5721520 0.4764767 0.5058857 0.5497430 
        VD23V     VD24N     VD24V     VD22N     VD22V     VD32N     VD32V     VD33N     VD33V 
    0.5609019 0.5333486 0.5137917 0.5043085 0.5135718 0.5188063 0.5048712 0.4710601 0.4498058 
        VD37N     VD37V     VD39N     VD39V     VD38N     VD38V     VD40N     VD40V     VD44N 
    0.4791925 0.4419878 0.4812183 0.5384752 0.5562247 0.5591057 0.5545795 0.5764205 0.5135818 
        VD44V     VD49N     VD49V 
    0.5532398 0.5006025 0.4919413 

    > data.norm<- neqc(data)
    > dim(data.norm)
    [1] 48803    30
    > expressed<-rowSums(data.norm$other$Detection_Pval < 0.05) >=3
    > data_expressed<-data.norm[expressed,]
    > dim(data_expressed)
    [1] 22249    30
    > tissue <-factor(targets$Tissue)
    > patient<-factor(targets$Sample_Name)
    > design <- model.matrix(~0+tissue+patient)
    > colnames(design)
     [1] "tissueL"     "tissueNL"    "patientVD18" "patientVD19" "patientVD20" "patientVD22"
     [7] "patientVD23" "patientVD24" "patientVD32" "patientVD33" "patientVD37" "patientVD38"
    [13] "patientVD39" "patientVD40" "patientVD44" "patientVD49"
    > cm<-makeContrasts(tissueL-tissueNL,levels = design)

    > cm
                 Contrasts
    Levels        tissue
      tissueL          1
      tissueNL        -1
      patientVD18      0
      patientVD19      0
      patientVD20      0
      .
      .  
      patientVD49      0

    > fit<-lmFit(data_expressed,design)
    > fit1<-contrasts.fit(fit, contrasts = cm)
    > fit2<-eBayes(fit)
    > tab <- topTable(fit2, coef=tissue)
    > View(tab)
    > tab <- topTable(fit2,adjust.method = "fdr", n= 100000)
    > colnames(tab)
     [1] "TargetID"    "tissueL"     "tissueNL"    "patientVD18" "patientVD19" "patientVD20"
     [7] "patientVD22" "patientVD23" "patientVD24" "patientVD32" "patientVD33" "patientVD37"
    [13] "patientVD38" "patientVD39" "patientVD40" "patientVD44" "patientVD49" "AveExpr"    
    [19] "F"           "P.Value"     "adj.P.Val"

[Benn]

So what is the error that you get? Please explain clearly what the issue is here.

[iti.gupta]

I am not getting any error But unfortunately, I am not getting the desired output. I want a single logFC and pvalue for all the 15 paired samples Also, I know some genes which for sure behave upregulated/downregulated as per the trend. Which is not the case here

[Benn]

Okay, I understand your issue now. Please try if this helps.

tab <- topTable(fit2, coef = 1, adjust.method = "fdr", n = 100000)

[russhh]

You've used fit in your call to eBayes, rather that the fit1 object you made while running contrasts.fit ... fit contains a coefficient for every coefficient in your design matrix, so running topTable over fit will do an F-Test over all the coefficients used in your model design. You don't want that.