limmaTestAndDesign.Rmdhttps://www.bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf
Method 1
#####################################
# Two groups
#####################################
design <- cbind(WT = 1, MUvsWT = c(0, 0, 1, 1, 1))
design## WT MUvsWT
## [1,] 1 0
## [2,] 1 0
## [3,] 1 1
## [4,] 1 1
## [5,] 1 1# Or by
Group <- factor(c(0, 0, 1, 1, 1))
design <- model.matrix(~Group)
colnames(design) <- c("WT", "MUvsWT")
design## WT MUvsWT
## 1 1 0
## 2 1 0
## 3 1 1
## 4 1 1
## 5 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$Group
## [1] "contr.treatment"Method 2
#####################################
# Two groups
#####################################
design <- cbind(WT = c(1, 1, 0, 0, 0), MU = c(0, 0, 1, 1, 1))
design## WT MU
## [1,] 1 0
## [2,] 1 0
## [3,] 0 1
## [4,] 0 1
## [5,] 0 1## WT MU
## 1 1 0
## 2 1 0
## 3 0 1
## 4 0 1
## 5 0 1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Group
## [1] "contr.treatment"## Contrasts
## Levels MUvsWT
## WT -1
## MU 1# fit <- lmFit(eset, design)
# fit2 <- contrasts.fit(fit, cont.matrix)
# fit2 <- eBayes(fit2)
# topTable(fit2, adjust="BH")
#####################################
# Three groups
#####################################
targets <- c("RNA1", "RNA1", "RNA2", "RNA2", "RNA3", "RNA3")
f <- factor(targets, levels = c("RNA1", "RNA2", "RNA3"))
design <- model.matrix(~ 0 + f)
colnames(design) <- c("RNA1", "RNA2", "RNA3")
design## RNA1 RNA2 RNA3
## 1 1 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 1 0
## 5 0 0 1
## 6 0 0 1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$f
## [1] "contr.treatment"contrast.matrix <- makeContrasts(RNA2 - RNA1, RNA3 - RNA2, RNA3 - RNA1, levels = design)
contrast.matrix## Contrasts
## Levels RNA2 - RNA1 RNA3 - RNA2 RNA3 - RNA1
## RNA1 -1 0 -1
## RNA2 1 -1 0
## RNA3 0 1 1sd <- 0.3 * sqrt(4 / rchisq(100, df = 4))
y <- matrix(rnorm(100 * 10, sd = sd), 100, 10)
rownames(y) <- paste("Gene", 1:100)
y[1:2, 4:6] <- y[1:2, 4:6] + 2# Method 1
Treat <- factor(c("C", "C", "C", "T", "T", "T", "G", "G", "G", "G"), levels = c("C", "T", "G"))
design <- model.matrix(~Treat)
design## (Intercept) TreatT TreatG
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 1 0
## 5 1 1 0
## 6 1 1 0
## 7 1 0 1
## 8 1 0 1
## 9 1 0 1
## 10 1 0 1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
TreatThreeT1 <- topTable(fit, coef = "TreatT")
TreatThreeG1 <- topTable(fit, coef = "TreatG")
contrast.matrix <- makeContrasts(TreatGT = TreatG - TreatT, levels = design)## Warning in makeContrasts(TreatGT = TreatG - TreatT, levels = design):
## Renaming (Intercept) to Intercept## Warning in contrasts.fit(fit, contrast.matrix): row names of contrasts
## don't match col names of coefficientsfit2 <- eBayes(fit2)
TreatThreeGT1 <- topTable(fit2, coef = "TreatGT")
# Method 2
design <- model.matrix(~ 0 + Treat)
design## TreatC TreatT TreatG
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 0 1 0
## 5 0 1 0
## 6 0 1 0
## 7 0 0 1
## 8 0 0 1
## 9 0 0 1
## 10 0 0 1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(TreatTC = TreatT - TreatC, TreatGC = TreatG - TreatC, TreatGT = TreatG - TreatT, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
TreatThreeT2 <- topTable(fit2, coef = "TreatTC")
TreatThreeG2 <- topTable(fit2, coef = "TreatGC")
TreatThreeGT2 <- topTable(fit2, coef = "TreatGT")
all.equal(TreatThreeT1, TreatThreeT2) # TRUE## [1] TRUE## [1] TRUE## [1] TRUEUsing design factor and contrast (1) A vs Other; (2) A vs (B+C)/2 are different. Method (2) is recommanded in most case as (1) used variation of all samples in Other group while (2) used the intra-group correlation for each individual group for the SEM.
See https://support.bioconductor.org/p/26251/ for more details
Treat <- factor(c("C", "C", "C", "T", "T", "T", "G", "G", "G", "G"), levels = c("C", "T", "G"))
design <- model.matrix(~Treat)
design## (Intercept) TreatT TreatG
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 1 0
## 5 1 1 0
## 6 1 1 0
## 7 1 0 1
## 8 1 0 1
## 9 1 0 1
## 10 1 0 1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
contrast.matrix <- makeContrasts(TreatGT = TreatG - TreatT, levels = design)## Warning in makeContrasts(TreatGT = TreatG - TreatT, levels = design):
## Renaming (Intercept) to Intercept## Warning in contrasts.fit(fit, contrast.matrix): row names of contrasts
## don't match col names of coefficientsfit2 <- eBayes(fit2)
TreatThreeGT1 <- topTable(fit2, coef = "TreatGT")
ySub=y[,-c(1)]
Treat <- factor(c("C", "C", "T", "T", "T", "G", "G", "G", "G"), levels = c("C", "T", "G"))
design <- model.matrix(~Treat)
design## (Intercept) TreatT TreatG
## 1 1 0 0
## 2 1 0 0
## 3 1 1 0
## 4 1 1 0
## 5 1 1 0
## 6 1 0 1
## 7 1 0 1
## 8 1 0 1
## 9 1 0 1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"fit <- lmFit(ySub, design)
fit <- eBayes(fit)
contrast.matrix <- makeContrasts(TreatGT = TreatG - TreatT, levels = design)## Warning in makeContrasts(TreatGT = TreatG - TreatT, levels = design):
## Renaming (Intercept) to Intercept## Warning in contrasts.fit(fit, contrast.matrix): row names of contrasts
## don't match col names of coefficientsfit2 <- eBayes(fit2)
TreatThreeGT2 <- topTable(fit2, coef = "TreatGT")
all.equal(TreatThreeGT1, TreatThreeGT2) # FALSE## [1] "Attributes: < Component \"row.names\": 8 string mismatches >"
## [2] "Component \"logFC\": Mean relative difference: 0.6315765"
## [3] "Component \"AveExpr\": Mean relative difference: 0.6534703"
## [4] "Component \"t\": Mean relative difference: 0.6070479"
## [5] "Component \"P.Value\": Mean relative difference: 0.1786546"
## [6] "Component \"adj.P.Val\": Mean relative difference: 0.1641983"
## [7] "Component \"B\": Mean relative difference: 0.07782985"# Method 1
Treat <- factor(c("C", "C", "C", "T", "T", "T", "C", "C", "T", "T"), levels = c("C", "T"))
Patient <- factor(c("A", "B", "C", "A", "B", "C", "D", "E", "D", "E"), levels = c("A", "B", "C", "D", "E"))
design <- model.matrix(~ Treat + Patient)
design## (Intercept) TreatT PatientB PatientC PatientD PatientE
## 1 1 0 0 0 0 0
## 2 1 0 1 0 0 0
## 3 1 0 0 1 0 0
## 4 1 1 0 0 0 0
## 5 1 1 1 0 0 0
## 6 1 1 0 1 0 0
## 7 1 0 0 0 1 0
## 8 1 0 0 0 0 1
## 9 1 1 0 0 1 0
## 10 1 1 0 0 0 1
## attr(,"assign")
## [1] 0 1 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"
##
## attr(,"contrasts")$Patient
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
TreatT1 <- topTable(fit, coef = "TreatT")
# Method 2
design <- model.matrix(~ 0 + Treat + Patient)
design## TreatC TreatT PatientB PatientC PatientD PatientE
## 1 1 0 0 0 0 0
## 2 1 0 1 0 0 0
## 3 1 0 0 1 0 0
## 4 0 1 0 0 0 0
## 5 0 1 1 0 0 0
## 6 0 1 0 1 0 0
## 7 1 0 0 0 1 0
## 8 1 0 0 0 0 1
## 9 0 1 0 0 1 0
## 10 0 1 0 0 0 1
## attr(,"assign")
## [1] 1 1 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"
##
## attr(,"contrasts")$Patient
## [1] "contr.treatment"fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(TreatTC = TreatT - TreatC, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
TreatT2 <- topTable(fit2, coef = "TreatTC")
all.equal(TreatT1, TreatT2) # TRUE## [1] TRUE########################################
# Two Factors
########################################
# Method 1
Treat <- factor(c("C", "C", "C", "T", "T", "T", "C", "C", "T", "T"), levels = c("C", "T"))
Mutation <- factor(c("Mut", "Mut", "Mut", "Mut", "Mut", "Mut", "WT", "WT", "WT", "WT"), levels = c("WT", "Mut"))
design <- model.matrix(~ Treat + Mutation)
design## (Intercept) TreatT MutationMut
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 1 1 1
## 5 1 1 1
## 6 1 1 1
## 7 1 0 0
## 8 1 0 0
## 9 1 1 0
## 10 1 1 0
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"
##
## attr(,"contrasts")$Mutation
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
MutationMut1 <- topTable(fit, coef = "MutationMut")
TreatTC1 <- topTable(fit, coef = "TreatT")
# Method 2
design <- model.matrix(~ 0 + Treat + Mutation)
design## TreatC TreatT MutationMut
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 0 1 1
## 5 0 1 1
## 6 0 1 1
## 7 1 0 0
## 8 1 0 0
## 9 0 1 0
## 10 0 1 0
## attr(,"assign")
## [1] 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Treat
## [1] "contr.treatment"
##
## attr(,"contrasts")$Mutation
## [1] "contr.treatment"fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(MutationMutWT = MutationMut, TreatTC = TreatT - TreatC, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
MutationMut2 <- topTable(fit2, coef = "MutationMutWT")
TreatTC2 <- topTable(fit2, coef = "TreatTC")
all.equal(TreatTC1, TreatTC2) # TRUE## [1] TRUE## [1] TRUE########################################
# Two Factors
########################################
# Method 1
Mutation1 <- factor(c("WT", "WT", "WT", "Mut", "Mut", "Mut", "WT", "WT", "WT", "WT"), levels = c("WT", "Mut"))
Mutation2 <- factor(c("Mut", "Mut", "Mut", "Mut", "Mut", "Mut", "WT", "WT", "WT", "WT"), levels = c("WT", "Mut"))
design <- model.matrix(~ Mutation1 + Mutation2)
design## (Intercept) Mutation1Mut Mutation2Mut
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 1 1 1
## 5 1 1 1
## 6 1 1 1
## 7 1 0 0
## 8 1 0 0
## 9 1 0 0
## 10 1 0 0
## attr(,"assign")
## [1] 0 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Mutation1
## [1] "contr.treatment"
##
## attr(,"contrasts")$Mutation2
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
Mutation1Mut1 <- topTable(fit, coef = "Mutation1Mut")
Mutation2Mut1 <- topTable(fit, coef = "Mutation2Mut")
# Method 2
design <- model.matrix(~ 0 + Mutation1 + Mutation2)
design## Mutation1WT Mutation1Mut Mutation2Mut
## 1 1 0 1
## 2 1 0 1
## 3 1 0 1
## 4 0 1 1
## 5 0 1 1
## 6 0 1 1
## 7 1 0 0
## 8 1 0 0
## 9 1 0 0
## 10 1 0 0
## attr(,"assign")
## [1] 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Mutation1
## [1] "contr.treatment"
##
## attr(,"contrasts")$Mutation2
## [1] "contr.treatment"fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(Mutation1Mut = Mutation1Mut - Mutation1WT, Mutation2Mut = Mutation2Mut, levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
Mutation1Mut2 <- topTable(fit2, coef = "Mutation1Mut")
Mutation2Mut2 <- topTable(fit2, coef = "Mutation2Mut")
all.equal(Mutation1Mut1, Mutation1Mut2) # TRUE## [1] TRUE## [1] TRUEMethod 2 is good and easy to write interactions
Similar as usersguide 9.5 Interaction Models: 2 × 2 Factorial Designs. But with data to make it easier to reproduce.
TS <- paste(Mutation, Treat, sep = ".")
TS <- factor(TS, levels = c("Mut.C", "Mut.T", "WT.C", "WT.T"))
design <- model.matrix(~ 0 + TS)
colnames(design) <- levels(TS)
design## Mut.C Mut.T WT.C WT.T
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 0 1 0 0
## 5 0 1 0 0
## 6 0 1 0 0
## 7 0 0 1 0
## 8 0 0 1 0
## 9 0 0 0 1
## 10 0 0 0 1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$TS
## [1] "contr.treatment"fit <- lmFit(y, design)
contrast.matrix <- makeContrasts(
TvsCinWT = WT.T - WT.C,
TvsCinMu = Mut.T - Mut.C,
MuvsWtinC = Mut.C - WT.C,
MuvsWtinT = Mut.T - WT.T,
Diff = (Mut.T - Mut.C) - (WT.T - WT.C),
levels = design
)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
TvsCinMu1 <- topTable(fit2, coef = "TvsCinMu")
TvsCinWT1 <- topTable(fit2, coef = "TvsCinWT")
MuvsWtinC1 <- topTable(fit2, coef = "MuvsWtinC")
MuvsWtinT1 <- topTable(fit2, coef = "MuvsWtinT")
diff1 <- topTable(fit2, coef = "Diff")# design <- model.matrix(~Mutation+Treat+Mutation:Treat)
design <- model.matrix(~ Mutation + Mutation:Treat)
design## (Intercept) MutationMut MutationWT:TreatT MutationMut:TreatT
## 1 1 1 0 0
## 2 1 1 0 0
## 3 1 1 0 0
## 4 1 1 0 1
## 5 1 1 0 1
## 6 1 1 0 1
## 7 1 0 0 0
## 8 1 0 0 0
## 9 1 0 1 0
## 10 1 0 1 0
## attr(,"assign")
## [1] 0 1 2 2
## attr(,"contrasts")
## attr(,"contrasts")$Mutation
## [1] "contr.treatment"
##
## attr(,"contrasts")$Treat
## [1] "contr.treatment"fit <- lmFit(y, design)
fit <- eBayes(fit)
TvsCinWT2 <- topTable(fit, coef = "MutationWT:TreatT")
TvsCinMu2 <- topTable(fit, coef = "MutationMut:TreatT")
all.equal(TvsCinWT1, TvsCinWT2) # TRUE## [1] TRUE## [1] TRUE# we could extract the interaction contrast Diff considered in "Single Factor" part by
fit2 <- contrasts.fit(fit, c(0, 0, -1, 1))
fit2 <- eBayes(fit2)
diff2 <- topTable(fit2)
all.equal(diff1, diff2) # TRUE## [1] TRUEdesign <- model.matrix(~ Mutation * Treat)
colnames(design) <- gsub(":|\\(|\\)", ".", colnames(design))
design## .Intercept. MutationMut TreatT MutationMut.TreatT
## 1 1 1 0 0
## 2 1 1 0 0
## 3 1 1 0 0
## 4 1 1 1 1
## 5 1 1 1 1
## 6 1 1 1 1
## 7 1 0 0 0
## 8 1 0 0 0
## 9 1 0 1 0
## 10 1 0 1 0
## attr(,"assign")
## [1] 0 1 2 3
## attr(,"contrasts")
## attr(,"contrasts")$Mutation
## [1] "contr.treatment"
##
## attr(,"contrasts")$Treat
## [1] "contr.treatment"# We need to use contrasts to extract all the comparisons of interest, as we are interested in Mu.S-Mu.U but it is not in design
fit <- lmFit(y, design)
# contrast.matrix <- cbind(TvsCinWT=c(0,0,1,0),TvsCinMu=c(0,0,1,1),Diff=c(0,0,0,1))
contrast.matrix <- makeContrasts(
MuvsWtinC = MutationMut,
MuvsWtinT = MutationMut + MutationMut.TreatT,
TvsCinWT = TreatT,
TvsCinMu = TreatT + MutationMut.TreatT,
Diff = MutationMut.TreatT,
levels = design
)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
MuvsWtinC3 <- topTable(fit2, coef = "MuvsWtinC")
MuvsWtinT3 <- topTable(fit2, coef = "MuvsWtinT")
TvsCinWT3 <- topTable(fit2, coef = "TvsCinWT")
TvsCinMu3 <- topTable(fit2, coef = "TvsCinMu")
diff3 <- topTable(fit2, coef = "Diff")
all.equal(TvsCinWT2, TvsCinWT3)## [1] TRUE## [1] TRUE## [1] TRUE## [1] TRUE## [1] TRUErawDataTable <- y
colnames(rawDataTable) <- paste0("Sample", 1:10)
makeDesignTable <- function(rawDataTable, factorTable,factorInDesign=NULL) {
if (is.null(factorInDesign)) {
factorInDesign=as.character(unique(factorTable[, 2]))
}
designTable <- matrix(0, nrow = ncol(rawDataTable), ncol = length(factorInDesign))
row.names(designTable) <- colnames(rawDataTable)
colnames(designTable) <- factorInDesign
designTable[as.matrix(factorTable)] <- 1
designTable
designTable <- cbind(Intercept = 1, designTable)
return(designTable)
}
performLimma <- function(rawDataTable, designTable, contrastsText = NULL) {
if (is.null(contrastsText)) {
contrastsText <- c(colnames(designTable)[-1])
}
fit <- lmFit(rawDataTable, designTable)
contrast.matrix <- makeContrasts(contrasts = contrastsText, levels = designTable)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
return(fit2)
}factorTable <- data.frame(Factor1 = c("Sample4", "Sample5", "Sample6", "Sample7", "Sample8", "Sample9", "Sample10"), Factor2 = c("TreatT", "TreatT", "TreatT", "TreatG", "TreatG", "TreatG", "TreatG"), stringsAsFactors = FALSE)
factorTable## Factor1 Factor2
## 1 Sample4 TreatT
## 2 Sample5 TreatT
## 3 Sample6 TreatT
## 4 Sample7 TreatG
## 5 Sample8 TreatG
## 6 Sample9 TreatG
## 7 Sample10 TreatGdesignTable <- makeDesignTable(rawDataTable, factorTable)
fit3 <- performLimma(rawDataTable, designTable)
TreatThreeT3 <- topTable(fit3, coef = "TreatT")
TreatThreeG3 <- topTable(fit3, coef = "TreatG")
all.equal(TreatThreeT1, TreatThreeT3) # True## [1] TRUE## [1] TRUEfactorTable <- data.frame(Factor1 = c("Sample4", "Sample5", "Sample6", "Sample9", "Sample10", "Sample2", "Sample3", "Sample5", "Sample6", "Sample7", "Sample8", "Sample9", "Sample10"), Factor2 = c("TreatT", "TreatT", "TreatT", "TreatT", "TreatT", "PatientB", "PatientC", "PatientB", "PatientC", "PatientD", "PatientE", "PatientD", "PatientE"), stringsAsFactors = FALSE)
factorTable## Factor1 Factor2
## 1 Sample4 TreatT
## 2 Sample5 TreatT
## 3 Sample6 TreatT
## 4 Sample9 TreatT
## 5 Sample10 TreatT
## 6 Sample2 PatientB
## 7 Sample3 PatientC
## 8 Sample5 PatientB
## 9 Sample6 PatientC
## 10 Sample7 PatientD
## 11 Sample8 PatientE
## 12 Sample9 PatientD
## 13 Sample10 PatientEdesignTable <- makeDesignTable(rawDataTable, factorTable)
fit3 <- performLimma(rawDataTable, designTable)
TreatT3 <- topTable(fit3, coef = "TreatT")
all.equal(TreatT2, TreatT3) # True## [1] TRUEfactorTable <- data.frame(Factor1 = c("Sample4", "Sample5", "Sample6", "Sample9", "Sample10", "Sample1", "Sample2", "Sample3", "Sample4", "Sample5", "Sample6"), Factor2 = c("TreatT", "TreatT", "TreatT", "TreatT", "TreatT", "MutationMut", "MutationMut", "MutationMut", "MutationMut", "MutationMut", "MutationMut"), stringsAsFactors = FALSE)
factorTable## Factor1 Factor2
## 1 Sample4 TreatT
## 2 Sample5 TreatT
## 3 Sample6 TreatT
## 4 Sample9 TreatT
## 5 Sample10 TreatT
## 6 Sample1 MutationMut
## 7 Sample2 MutationMut
## 8 Sample3 MutationMut
## 9 Sample4 MutationMut
## 10 Sample5 MutationMut
## 11 Sample6 MutationMutdesignTable <- makeDesignTable(rawDataTable, factorTable)
fit3 <- performLimma(rawDataTable, designTable)
TreatTC3 <- topTable(fit3, coef = "TreatT")
MutationMut3 <- topTable(fit3, coef = "MutationMut")
all.equal(TreatTC2, TreatTC3) # True## [1] TRUE## [1] TRUE########################################################
# Data prepare. For test, Should NOT in pipeline
########################################################
# library(reshape2)
library(tidyverse)
setwd("D:/OneDriveWork/OneDriveVanderbilt/work/AnthonyDaniels/20190509_proteomicsReDo")
proteinExpression <- readRDS("20190827_withMatchedNormalProteinExpression.rds")
# rawData=proteinExpression %>% spread(Sample, nLog2)
# order by abs value of rev(nLog2) so that distinct will keep the first (largest) row
temp <- abs(proteinExpression$nLog2)
rawData <- proteinExpression[rev(order(temp)), ] %>%
distinct(pID, Sample, .keep_all = TRUE) %>%
spread(Sample, nLog2)
rawData <- rawData %>% distinct(pID, .keep_all = TRUE) # remove duplicate pID, most are error
row.names(rawData) <- rawData[, 1]
rawData <- rawData[, -1]
write.csv(rawData, "D:\\OneDriveWork\\OneDriveVanderbilt\\work\\AnthonyDaniels\\20190903_proteomicsLimma\\20190903_threeTumorNormalPair.csv")