batchEffectAnalysis.rmd

```{r data, echo = FALSE}
inpD <-"F:\\promec\\Animesh\\HUNT\\txt106dpMBR\\"
data <- read.delim(paste0(inpD,"proteinGroups.txt"),row.names=1,sep="\t",header = T)
summary(data)
decoyPrefix="REV"
contaminantPrefix="CON"
dataC="LFQ.intensity."
```


```{r dataNorm, echo = FALSE}
dataNorm=log2(data[,grep(dataC, names(data))])
#set.seed(1)
#dataNorm[dataNorm==0]<-NaN#rnorm(1,mean=mean(dataNorm),sd=sd(dataNorm))
dataNorm[dataNorm==-Inf]=NA
summary(dataNorm)
```

```{r write-output, echo = FALSE}
write.table(dataNorm,file=paste0(inpD,"log2data.txt"), sep = "\t")
#dump(dataNorm,file=paste0(inpD,"dataNorm.R"))
```


```{r impute, echo = FALSE}
#install.packages('mice')
library(mice)
#install.packages('randomForest')
library(randomForest)
dataNormImp=mice(dataNorm, method="rf")
dataNormImpCom <- complete(dataNormImp,1)
row.names(dataNormImpCom)<-row.names(dataNorm)
summary(dataNormImpCom)
```

```{r write-output, echo = FALSE}
#write.csv(dataNormImpCom,file=paste0(inpD,"log2dataImp.csv"))
#write.csv(factors,file=paste0(inpD,"dataNormImpComFactor.csv"))
dataNormImpCom <- read.csv(paste0(inpD,"log2dataImp.csv"),row.names=1,header = T)
#factors<-read.csv(paste0(inpD,"dataNormImpComFactor.csv"))
#dump(dataNorm,file=paste0(inpD,"dataNorm.R"))
```

```{r PCA, echo = FALSE}
plot(princomp(dataNormImpCom))
biplot(prcomp(dataNormImpCom,scale = F))
biplot(princomp(as.matrix(t(dataNormImpCom)),scale = T))
```


```{r proBatch}
#https://github.com/symbioticMe/proBatch
#needs R version 3.6 
#install.packages("BiocManager")
#BiocManager::install("proBatch")
#install.packages("devtools")
library(devtools)
#install_github("symbioticMe/proBatch", build_vignettes = TRUE)
library(proBatch)
help(proBatch)
```


```{r condition, echo = FALSE}
factors=as.factor(sapply(strsplit(names(data)[grep(dataC, names(data))],"[0-9]+"),'[[', 2))
factors=as.integer(factors)
#factors=as.character(sapply(strsplit(names(data)[grep(dataC, names(data))],"[0-9]+"),'[[', 2))
```

```{r DEqMS}
#https://rdrr.io/bioc/DEqMS/f/vignettes/DEqMS-package-vignette.Rmd
#install.packages('devtools')
#devtools::install_github("syspremed/exploBATCH",verbose=TRUE)
#BiocManager::install("DEqMS")
library(DEqMS)
#url <- "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2016/06/PXD004163/Yan_miR_Protein_table.flatprottable.txt"
#download.file(url, destfile = "./miR_Proteintable.txt",method = "auto")
df.prot = read.table("F:\\promec\\Animesh\\Yan_miR_Protein_table.flatprottable.txt",stringsAsFactors = FALSE,header = TRUE, quote = "", comment.char = "",sep = "\t")
# filter at 1% protein FDR and extract TMT quantifications
TMT_columns = seq(15,33,2)
dat = df.prot[df.prot$miR.FASP_q.value<0.01,TMT_columns]
rownames(dat) = df.prot[df.prot$miR.FASP_q.value<0.01,]$Protein.accession
# The protein dataframe is a typical protein expression matrix structure
# Samples are in columns, proteins are in rows
# use unique protein IDs for rownames
# to view the whole data frame, use the command View(dat)
dat.log = log2(dat)
#remove rows with NAs
dat.log = na.omit(dat.log)
boxplot(dat.log,las=2,main="TMT10plex data PXD004163")
# if there is only one factor, such as treatment. You can define a vector with
# the treatment group in the same order as samples in the protein table.
cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl",
"miR372","miR519","ctrl","miR191","miR372"))
# The function model.matrix is used to generate the design matrix
design = model.matrix(~0+cond) # 0 means no intercept for the linear model
colnames(design) = gsub("cond","",colnames(design))
# you can define one or multiple contrasts here
x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl",
       "miR372-miR519","miR372-miR191","miR519-miR191")
contrast =  makeContrasts(contrasts=x,levels=design)
fit1 <- lmFit(dat.log, design)
fit2 <- contrasts.fit(fit1,contrasts = contrast)
fit3 <- eBayes(fit2)
# assign a extra variable `count` to fit3 object, telling how many PSMs are 
# quantifed for each protein
library(matrixStats)
count_columns = seq(16,34,2)
psm.count.table = data.frame(count = rowMins(
  as.matrix(df.prot[,count_columns])), row.names =  df.prot$Protein.accession)
fit3$count = psm.count.table[rownames(fit3$coefficients),"count"]
fit4 = spectraCounteBayes(fit3)
# n=30 limits the boxplot to show only proteins quantified by <= 30 PSMs.
VarianceBoxplot(fit4,n=30,main="TMT10plex dataset PXD004163",xlab="PSM count")
VarianceScatterplot(fit4,main="TMT10plex dataset PXD004163")
DEqMS.results = outputResult(fit4,coef_col = 1)
#if you are not sure which coef_col refers to the specific contrast,type
head(fit4$coefficients)
# a quick look on the DEqMS results table
head(DEqMS.results)
# Save it into a tabular text file
write.table(DEqMS.results,paste0(inpD,"DEqMS.results.miR372-ctrl.txt"),sep = "\t",row.names = F,quote=F)
install.packages("ggrepel")
library(ggrepel)
# Use ggplot2 allows more flexibility in plotting
DEqMS.results$log.sca.pval = -log10(DEqMS.results$sca.P.Value)
ggplot(DEqMS.results, aes(x = logFC, y =log.sca.pval )) + 
    geom_point(size=0.5 )+
    theme_bw(base_size = 16) + # change theme
    xlab(expression("log2(miR372/ctrl)")) + # x-axis label
    ylab(expression(" -log10(P-value)")) + # y-axis label
    geom_vline(xintercept = c(-1,1), colour = "red") + # Add fold change cutoffs
    geom_hline(yintercept = 3, colour = "red") + # Add significance cutoffs
    geom_vline(xintercept = 0, colour = "black") + # Add 0 lines
    scale_colour_gradient(low = "black", high = "black", guide = FALSE)+
    geom_text_repel(data=subset(DEqMS.results, abs(logFC)>1&log.sca.pval > 3),
                    aes( logFC, log.sca.pval ,label=gene)) # add gene label

fit4$p.value = fit4$sca.p
# volcanoplot highlight top 20 proteins ranked by p-value here
volcanoplot(fit4,coef=1, style = "p-value", highlight = 20,names=rownames(fit4$coefficients))


design = model.matrix(~0+as.factor(factors))
colnames(design) = gsub("as\\.factor\\(factors\\)","",colnames(design))
contrast =  makeContrasts(contrasts=colnames(dataNormImpCom),levels=design)
contrast =  makeContrasts(contrasts=colnames(design),levels=design)
fit1 <- lmFit(dat.log, design)
fit2 <- contrasts.fit(fit1,contrasts = contrast)
fit3 <- eBayes(fit2)

df.LFQ<-dataNormImpCom
boxplot(df.LFQ)
library(matrixStats)
pep.count.table = data.frame(count = rowMins(as.matrix(data[,112:217])))
# Minimum peptide count of some proteins can be 0
# add pseudocount 1 to all proteins
pep.count.table$count = pep.count.table$count+1
class<-as.factor(factors)
design = model.matrix(~0+class) 

protein.matrix = as.matrix(df.LFQ)
fit1 = lmFit(protein.matrix,design = design)
?makeContrasts
cont <- makeContrasts(c(1)-c(2),levels = design)
cont <- makeContrasts(c(1),levels = design)
fit2 = contrasts.fit(fit1,contrasts = cont)
fit3 <- eBayes(fit2)

fit3$count = pep.count.table$count#[rownames(fit3$coefficients),"count"]

#check the values in the vector fit3$count
#if min(fit3$count) return NA or 0, you should troubleshoot the error first
min(fit3$count)

fit4 = spectraCounteBayes(fit3)
VarianceBoxplot(fit4, n=20, main = "lfqHUNT",xlab="peptide count + 1")
```

```{r exploBatch, echo = FALSE}
Download and Read the input protein table
url <- "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2016/06/PXD004163/Yan_miR_Protein_table.flatprottable.txt"
download.file(url, destfile = "./miR_Proteintable.txt",method = "auto")

df.prot = read.table("miR_Proteintable.txt",stringsAsFactors = FALSE,
                     header = TRUE, quote = "", comment.char = "",sep = "\t")
Extract quant data columns for DEqMS
# filter at 1% protein FDR and extract TMT quantifications
TMT_columns = seq(15,33,2)
dat = df.prot[df.prot$miR.FASP_q.value<0.01,TMT_columns]
rownames(dat) = df.prot[df.prot$miR.FASP_q.value<0.01,]$Protein.accession
# The protein dataframe is a typical protein expression matrix structure
# Samples are in columns, proteins are in rows
# use unique protein IDs for rownames
# to view the whole data frame, use the command View(dat)
If the protein table is relative abundance (ratios) or intensity values, Log2 transform the data. Systematic effects and variance components are usually assumed to be additive on log scale (Oberg AL. et al JPR 2008; Hill EG. et al JPR 2008).

dat.log = log2(dat)
#remove rows with NAs
dat.log = na.omit(dat.log)
Use boxplot to check if the samples have medians centered. if not, do median centering.

boxplot(dat.log,las=2,main="TMT10plex data PXD004163")
# Here the data is already median centered, we skip the following step. 
# dat.log = equalMedianNormalization(dat.log)
Make design table.
A design table is used to tell how samples are arranged in different groups/classes.

# if there is only one factor, such as treatment. You can define a vector with
# the treatment group in the same order as samples in the protein table.
cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl",
"miR372","miR519","ctrl","miR191","miR372"))

# The function model.matrix is used to generate the design matrix
design = model.matrix(~0+cond) # 0 means no intercept for the linear model
colnames(design) = gsub("cond","",colnames(design))
Make contrasts
In addition to the design, you need to define the contrast, which tells the model to compare the differences between specific groups. Start with the Limma part.

# you can define one or multiple contrasts here
x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl",
       "miR372-miR519","miR372-miR191","miR519-miR191")
contrast =  makeContrasts(contrasts=x,levels=design)
fit1 <- lmFit(dat.log, design)
fit2 <- contrasts.fit(fit1,contrasts = contrast)
fit3 <- eBayes(fit2)
DEqMS analysis
The above shows Limma part, now we use the function spectraCounteBayes in DEqMS to correct bias of variance estimate based on minimum number of psms per protein used for quantification.We use the minimum number of PSMs used for quantification within and across experiments to model the relation between variance and PSM count.(See original paper)

# assign a extra variable `count` to fit3 object, telling how many PSMs are 
# quantifed for each protein
library(matrixStats)
count_columns = seq(16,34,2)
psm.count.table = data.frame(count = rowMins(
  as.matrix(df.prot[,count_columns])), row.names =  df.prot$Protein.accession)
fit3$count = psm.count.table[rownames(fit3$coefficients),"count"]
fit4 = spectraCounteBayes(fit3)
Outputs of spectraCounteBayes:
object is augmented form of "fit" object from eBayes in Limma, with the additions being:
sca.t - Spectra Count Adjusted posterior t-value
sca.p - Spectra Count Adjusted posterior p-value
sca.dfprior - DEqMS estimated prior degrees of freedom
sca.priorvar- DEqMS estimated prior variance
sca.postvar - DEqMS estimated posterior variance
model - fitted model

Visualize the fit curve - variance dependence on quantified PSM
# n=30 limits the boxplot to show only proteins quantified by <= 30 PSMs.
VarianceBoxplot(fit4,n=30,main="TMT10plex dataset PXD004163",xlab="PSM count")
VarianceScatterplot(fit4,main="TMT10plex dataset PXD004163")
Extract the results as a data frame and save it
DEqMS.results = outputResult(fit4,coef_col = 1)
#if you are not sure which coef_col refers to the specific contrast,type
head(fit4$coefficients)

# a quick look on the DEqMS results table
head(DEqMS.results)
# Save it into a tabular text file
write.table(DEqMS.results,"DEqMS.results.miR372-ctrl.txt",sep = "\t",
            row.names = F,quote=F)
Explaination of the columns in DEqMS.results:
logFC - log2 fold change between two groups, Here it's log2(miR372/ctrl).
AveExpr - the mean of the log2 ratios/intensities across all samples. Since input matrix is log2 ratio values, it is the mean log2 ratios of all samples.
t - Limma output t-statistics
P.Value- Limma p-values
adj.P.Val - BH method adjusted Limma p-values
B - Limma B values
count - PSM/peptide count values you assigned
sca.t - DEqMS t-statistics
sca.P.Value - DEqMS p-values
sca.adj.pval - BH method adjusted DEqMS p-values

Make volcanoplot
We recommend to plot p-values on y-axis instead of adjusted pvalue or FDR.
Read about why here.

library(ggrepel)
# Use ggplot2 allows more flexibility in plotting

DEqMS.results$log.sca.pval = -log10(DEqMS.results$sca.P.Value)
ggplot(DEqMS.results, aes(x = logFC, y =log.sca.pval )) + 
    geom_point(size=0.5 )+
    theme_bw(base_size = 16) + # change theme
    xlab(expression("log2(miR372/ctrl)")) + # x-axis label
    ylab(expression(" -log10(P-value)")) + # y-axis label
    geom_vline(xintercept = c(-1,1), colour = "red") + # Add fold change cutoffs
    geom_hline(yintercept = 3, colour = "red") + # Add significance cutoffs
    geom_vline(xintercept = 0, colour = "black") + # Add 0 lines
    scale_colour_gradient(low = "black", high = "black", guide = FALSE)+
    geom_text_repel(data=subset(DEqMS.results, abs(logFC)>1&log.sca.pval > 3),
                    aes( logFC, log.sca.pval ,label=gene)) # add gene label
you can also use volcanoplot function from Limma. However, it uses p.value from Limma. If you want to plot sca.pvalue from DEqMS, you need to modify the fit4 object.

fit4$p.value = fit4$sca.p
# volcanoplot highlight top 20 proteins ranked by p-value here
volcanoplot(fit4,coef=1, style = "p-value", highlight = 20,
            names=rownames(fit4$coefficients))
DEqMS analysis using MaxQuant outputs (label-free data)
Here we analyze a published label-free benchmark dataset in which either 10 or 30 µg of E. coli protein extract was spiked into human protein extracts (50 µg) in triplicates (Cox J et al MCP 2014). The data was searched by MaxQuant software and the output file "proteinGroups.txt" was used here.

url2 <- "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2014/09/PXD000279/proteomebenchmark.zip"
download.file(url2, destfile = "./PXD000279.zip",method = "auto")
unzip("PXD000279.zip")
Read protein table as input and filter it
df.prot = read.table("proteinGroups.txt",header=T,sep="\t",stringsAsFactors = F,
                        comment.char = "",quote ="")

# remove decoy matches and matches to contaminant
df.prot = df.prot[!df.prot$Reverse=="+",]
df.prot = df.prot[!df.prot$Contaminant=="+",]

# Extract columns of LFQ intensites
df.LFQ = df.prot[,89:94]
df.LFQ[df.LFQ==0] <- NA

rownames(df.LFQ) = df.prot$Majority.protein.IDs
df.LFQ$na_count_H = apply(df.LFQ,1,function(x) sum(is.na(x[1:3])))
df.LFQ$na_count_L = apply(df.LFQ,1,function(x) sum(is.na(x[4:6])))
# Filter protein table. DEqMS require minimum two values for each group.
df.LFQ.filter = df.LFQ[df.LFQ$na_count_H<2 & df.LFQ$na_count_L<2,1:6]
Make a data frame of unique peptide count per protein
library(matrixStats)
# we use minimum peptide count among six samples
# count unique+razor peptides used for quantification
pep.count.table = data.frame(count = rowMins(as.matrix(df.prot[,19:24])),
                             row.names = df.prot$Majority.protein.IDs)
# Minimum peptide count of some proteins can be 0
# add pseudocount 1 to all proteins
pep.count.table$count = pep.count.table$count+1
DEqMS analysis on LFQ data
protein.matrix = log2(as.matrix(df.LFQ.filter))

class = as.factor(c("H","H","H","L","L","L"))
design = model.matrix(~0+class) # fitting without intercept

fit1 = lmFit(protein.matrix,design = design)
cont <- makeContrasts(classH-classL, levels = design)
fit2 = contrasts.fit(fit1,contrasts = cont)
fit3 <- eBayes(fit2)

fit3$count = pep.count.table[rownames(fit3$coefficients),"count"]

#check the values in the vector fit3$count
#if min(fit3$count) return NA or 0, you should troubleshoot the error first
min(fit3$count)

fit4 = spectraCounteBayes(fit3)
Visualize the fit curve
VarianceBoxplot(fit4, n=20, main = "Label-free dataset PXD000279",
                xlab="peptide count + 1")
Extract outputs from DEqMS
DEqMS.results = outputResult(fit4,coef_col = 1)
# Add Gene names to the data frame
rownames(df.prot) = df.prot$Majority.protein.IDs
DEqMS.results$Gene.name = df.prot[DEqMS.results$gene,]$Gene.names
head(DEqMS.results)
write.table(DEqMS.results,"H-L.DEqMS.results.txt",sep = "\t",
            row.names = F,quote=F)
DEqMS analysis using a PSM table (isobaric labelled data)
If you want to try different methods to estimate protein abundance,you can start with a PSM table and use provided functions in DEqMS to summarize PSM quant data into protein quant data. Four different functions are included: medianSweeping,medianSummary,medpolishSummary,farmsSummary. Check PDF reference manual for detailed description.

Read PSM table input
### retrieve example PSM dataset from ExperimentHub
library(ExperimentHub)
eh = ExperimentHub()
query(eh, "DEqMS")
dat.psm = eh[["EH1663"]]
dat.psm.log = dat.psm
dat.psm.log[,3:12] =  log2(dat.psm[,3:12])
head(dat.psm.log)
Summarization and Normalization
Here, median sweeping is used to summarize PSMs intensities to protein log2 ratios. In this procedure, we substract the spectrum log2 intensity from the median log2 intensities of all samples. The relative abundance estimate for each protein is calculated as the median over all PSMs belonging to this protein.(Herbrich et al JPR 2012 and D'Angelo et al JPR 2016).
Assume the log2 intensity of PSM i in sample j is $y_{i,j}$, its relative log2 intensity of PSM i in sample j is $y'{i,j}$: $$y'{i,j} = y_{i,j} - median_{j'\in ctrl}\ y_{i,j'} $$ Relative abundance of protein k in sample j $Y_{k,j}$ is calculated as: $$Y_{k,j} = median_{i\in protein\ k}\ y'_{i,j} $$

Correction for differences in amounts of material loaded in the channels is then done by subtracting the channel median from the relative abundance (log2 ratio), centering all channels to have median log2 value of zero.

dat.gene.nm = medianSweeping(dat.psm.log,group_col = 2)
boxplot(dat.gene.nm,las=2,ylab="log2 ratio",main="TMT10plex dataset PXD004163")
DEqMS analysis
gene.matrix = as.matrix(dat.gene.nm)

# make design table
cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl",
"miR372","miR519","ctrl","miR191","miR372"))
design = model.matrix(~0+cond) 
colnames(design) = gsub("cond","",colnames(design))

#limma part analysis
fit1 <- lmFit(gene.matrix,design)
x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl")
contrast =  makeContrasts(contrasts=x,levels=design)
fit2 <- eBayes(contrasts.fit(fit1,contrasts = contrast))

#DEqMS part analysis
psm.count.table = as.data.frame(table(dat.psm$gene))
rownames(psm.count.table) = psm.count.table$Var1

fit2$count = psm.count.table[rownames(fit2$coefficients),2]
fit3 = spectraCounteBayes(fit2)
# extract DEqMS results
DEqMS.results = outputResult(fit3,coef_col = 1) 
head(DEqMS.results)
write.table(DEqMS.results,"DEqMS.results.miR372-ctrl.fromPSMtable.txt",
            sep = "\t",row.names = F,quote=F)
Generate variance ~ PMS count boxplot, check if the DEqMS correctly find the relation between prior variance and PSM count

VarianceBoxplot(fit3,n=20, xlab="PSM count",main="TMT10plex dataset PXD004163")
PSM/Peptide profile plot
Only possible if you read a PSM or peptide table as input. peptideProfilePlot function will plot log2 intensity of each PSM/peptide of the protein in the input table.

peptideProfilePlot(dat=dat.psm.log,col=2,gene="TGFBR2")
# col=2 is tell in which column of dat.psm.log to look for the gene
Comparing DEqMS to other methods
The following steps are not required for get the results from DEqMS. it is used to help users to understand the method better and the differences to other methods. Here we use the TMT labelled data PXD004163 as an example.

Compare the variance estimate in DEqMS and Limma
Prior variance comparison between DEqMS and Limma
VarianceScatterplot(fit3, xlab="log2(PSM count)")
limma.prior = fit3$s2.prior
abline(h = log(limma.prior),col="green",lwd=3 )
legend("topright",legend=c("DEqMS prior variance","Limma prior variance"),
        col=c("red","green"),lwd=3)
Residual plot for DEqMS and Limma
op <- par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0))
Residualplot(fit3,  xlab="log2(PSM count)",main="DEqMS")
x = fit3$count
y = log(limma.prior) - log(fit3$sigma^2)
plot(log2(x),y,ylim=c(-6,2),ylab="Variance(estimated-observed)", pch=20, cex=0.5,
     xlab = "log2(PSMcount)",main="Limma")
Posterior variance comparison between DEqMS and Limma
The plot here shows posterior variance of proteins "shrink" toward the fitted value to different extent depending on PSM number.

library(LSD)
op <- par(mfrow=c(1,2), mar=c(4,4,4,1), oma=c(0.5,0.5,0.5,0))
x = fit3$count
y = fit3$s2.post
heatscatter(log2(x),log(y),pch=20, xlab = "log2(PSMcount)", 
     ylab="log(Variance)",
     main="Posterior Variance in Limma")

y = fit3$sca.postvar
heatscatter(log2(x),log(y),pch=20, xlab = "log2(PSMcount)",
     ylab="log(Variance)", 
     main="Posterior Variance in DEqMS")
Compare p-values from DEqMS to ordinary t-test, ANOVA and Limma
We first apply t.test to detect significant protein changes between ctrl samples and miR372 treated samples, both have three replicates.

T-test analysis
pval.372 = apply(dat.gene.nm, 1, function(x) 
t.test(as.numeric(x[c(1,5,8)]), as.numeric(x[c(3,6,10)]))$p.value)

logFC.372 = rowMeans(dat.gene.nm[,c(3,6,10)])-rowMeans(dat.gene.nm[,c(1,5,8)])
Generate a data.frame of t.test results, add PSM count values and order the table by p-value.

ttest.results = data.frame(gene=rownames(dat.gene.nm),
                    logFC=logFC.372,P.Value = pval.372, 
                    adj.pval = p.adjust(pval.372,method = "BH")) 

ttest.results$PSMcount = psm.count.table[ttest.results$gene,"count"]
ttest.results = ttest.results[with(ttest.results, order(P.Value)), ]
head(ttest.results)
Anova analysis
Anova analysis is equivalent to linear model analysis. The difference to Limma analysis is that estimated variance is not moderated using empirical bayesian approach as it is done in Limma.

ord.t = fit1$coefficients[, 1]/fit1$sigma/fit1$stdev.unscaled[, 1]
ord.p = 2*pt(abs(ord.t), fit1$df.residual, lower.tail = FALSE)
ord.q = p.adjust(ord.p,method = "BH")
anova.results = data.frame(gene=names(fit1$sigma),
                            logFC=fit1$coefficients[,1],
                            t=ord.t, 
                            P.Value=ord.p, 
                            adj.P.Val = ord.q)

anova.results$PSMcount = psm.count.table[anova.results$gene,"count"]
anova.results = anova.results[with(anova.results,order(P.Value)),]

head(anova.results)
Limma
Extract limma results using topTable function, coef = 1 allows you to extract the specific contrast (miR372-ctrl), option n= Inf output all rows.

limma.results = topTable(fit2,coef = 1,n= Inf)
limma.results$gene = rownames(limma.results)
#Add PSM count values in the data frame
limma.results$PSMcount = psm.count.table[limma.results$gene,"count"]

head(limma.results)
Visualize the distribution of p-values by different analysis
plotting all proteins ranked by p-values.

plot(sort(-log10(limma.results$P.Value),decreasing = TRUE), 
    type="l",lty=2,lwd=2, ylab="-log10(p-value)",ylim = c(0,10),
    xlab="Proteins ranked by p-values",
    col="purple")
lines(sort(-log10(DEqMS.results$sca.P.Value),decreasing = TRUE), 
        lty=1,lwd=2,col="red")
lines(sort(-log10(anova.results$P.Value),decreasing = TRUE), 
        lty=2,lwd=2,col="blue")
lines(sort(-log10(ttest.results$P.Value),decreasing = TRUE), 
        lty=2,lwd=2,col="orange")
legend("topright",legend = c("Limma","DEqMS","Anova","t.test"),
        col = c("purple","red","blue","orange"),lty=c(2,1,2,2),lwd=2)
plotting top 500 proteins ranked by p-values.

plot(sort(-log10(limma.results$P.Value),decreasing = TRUE)[1:500], 
    type="l",lty=2,lwd=2, ylab="-log10(p-value)", ylim = c(2,10),
    xlab="Proteins ranked by p-values",
    col="purple")
lines(sort(-log10(DEqMS.results$sca.P.Value),decreasing = TRUE)[1:500], 
        lty=1,lwd=2,col="red")
lines(sort(-log10(anova.results$P.Value),decreasing = TRUE)[1:500], 
        lty=2,lwd=2,col="blue")
lines(sort(-log10(ttest.results$P.Value),decreasing = TRUE)[1:500], 
        lty=2,lwd=2,col="orange")
legend("topright",legend = c("Limma","DEqMS","Anova","t.test"),
        col = c("purple","red","blue","orange"),lty=c(2,1,2,2),lwd=2)


#install.packages('devtools')
#devtools::install_github("syspremed/exploBATCH",verbose=TRUE)
#require(exploBATCHbreast)
library("exploBATCH")
#?expBATCH
#data(Breast)
#data(batchBreast)
expBATCH(D=t(dataNormImpCom),batchCL=(factors),mindim=1,maxdim=7,method="ppcca")#,SDselect=2)
#expBATCH(D=Breast,batchCL=batchBreast,Conf=NA,mindim=5,maxdim=7,method="ppcca")#,SDselect=2)
#expBATCH(Breast,batchBreast)
#expBATCH(t(dataNormImpCom),factors)
```


```{r ROTS, echo = FALSE}
#install.packages("BiocManager")
#BiocManager::install("ROTS", version = "3.8")
dataNormImpCom<-read.delim("L:\\promec\\Elite\\LARS\\2014\\desember\\christiano\\plantaris-walker-hit.txt",row.names=1,sep="\t",header = T)
dataNormImpCom[is.na(dataNormImpCom)]=5
summary(dataNormImpCom)
factors<-c(rep(0,5), rep(1,5))
library(ROTS)
results = ROTS(data = dataNormImpCom, groups = factors , B = 1000 , K = 500 , seed = 1234)
write.csv(summary(results, fdr = 1),file=paste0(inpD,"L:\\promec\\Elite\\LARS\\2014\\desember\\christiano\\plantaris-walker-hit.rots.csv"), sep = "\t")

names(results) 
summary(results, fdr = 1)
plot(results, fdr = 0.5, type = "volcano")
#P61203-2;P61203
#plot(results, fdr = 0.5, type = "heatmap")
```



```{r MSstats, echo = FALSE}
#https://bioconductor.org/packages/3.3/bioc/vignettes/MSstats/inst/doc/MSstats-manual.pdf
library('MSstats')
QuantData<-dataProcess(SRMRawData)
head(QuantData$ProcessedData)
## based on multiple comparisons (T1 vs T3; T1 vs T7; T1 vs T9)
comparison1<-matrix(c(-1,0,1,0,0,0,0,0,0,0),nrow=1)
comparison2<-matrix(c(-1,0,0,0,0,0,1,0,0,0),nrow=1)
comparison3<-matrix(c(-1,0,0,0,0,0,0,0,1,0),nrow=1)
comparison<-rbind(comparison1,comparison2, comparison3)
row.names(comparison)<-c("T3-T1","T7-T1","T9-T1")
testResultMultiComparisons<-groupComparison(contrast.matrix=comparison,data=QuantData)
## Calculate sample size for future experiments:
#(1) Minimal number of biological replicates per condition
designSampleSize(data=testResultMultiComparisons$fittedmodel,numSample=TRUE,
desiredFC=c(1.25,1.75),FDR=0.05,power=0.8)
#(2) Power calculation
designSampleSize(data=testResultMultiComparisons$fittedmodel,numSample=2,
desiredFC=c(1.25,1.75),FDR=0.05,power=TRUE)```
```

```{r SAM, echo = FALSE}
#https://github.com/MikeJSeo/SAM
install.packages(c("samr", "matrixStats", "GSA", "shiny", "shinyFiles", "openxlsx"))
source("http://bioconductor.org/biocLite.R")
biocLite("impute")
library(shiny)
library(shinyFiles)
runGitHub("SAM", "MikeJSeo")
```




```{r NA}
#impute(data, method = "mixed",randna = fData(naset)$randna,mar = "knn", mnar = "min")
md.pattern(datasel)
dataimp <- mice(datasel,m=length(colnames(datasel))/2,maxit=length(colnames(datasel)),meth='pmm',seed=1)
dataimp <- complete(dataimp,length(colnames(datasel))/2)
heatmap(as.matrix(dataimp),col=colscl)
```

```{r Scale}
datascale=t(scale(t(dataimp))) # takes away the polynomial in YFP data
heatmap(as.matrix(datascale),col=colscl)
```

```{r DistCor}
dataasmat<-data.matrix(datascale)
cor(data)
dissimilarityc <- 1 - cor(dataasmat,method = "s")
distancec <- as.dist(dissimilarityc)
col.clus <- hclust(distancec, "aver")
dissimilarityr <- 1 - cor(t(dataasmat),method = "s")
distancer <- as.dist(dissimilarityr)
row.clus <- hclust(distancer, "aver")
heatmap(as.matrix(datascale),col=colscl)
cr <- rainbow(nrow(dataasmat))
cc <- rainbow(ncol(dataasmat))
heatmap((na.omit(dataasmat)), RowSideColors = cr, ColSideColors = cc,Rowv = as.dendrogram(row.clus), Colv = as.dendrogram(col.clus),cexRow=1.5,cexCol=1.5 )
```


```{r Quantile}
dataasmat<-normalize.quantiles(dataasmat)
rn<-rownames(dataasmat)
cn<-colnames(dataasmat)
colnames(dataasmat)<-cn
rownames(dataasmat)<-rn
```

```{r AdjPval}
p.adjust(10^(-data$X.log10P),method = c("hochberg"), n = length(data$X.log10P))
p.adjust(data$pv,method = c("hochberg"), n = length(data$pv))
```

```{r interpolate}
data<-read.table("c:/Users/animeshs/Downloads/Sunshine hours.txt",sep='\t',header=T)
summary(data$SONNENSCHEINDAUER)
dp=approx(data$SONNENSCHEINDAUER)
summary(dp)
plot(dp$SONNENSCHEINDAUER,data$SONNENSCHEINDAUER)
install.packages('zoo')
library('zoo')
dp<-rollmean(data,1)
dp<-na.approx(data)
summary(dp)
write.table(dp,"c:/Users/animeshs/Downloads/SunshineDP.txt",sep='\t')
```



```{r lmfit}
source("https://bioconductor.org/biocLite.R")
biocLite("golubEsets")
library(golubEsets)
data()
data(golubTrain)
golubTrain
golub<-data.frame(golubTrain@exprs)
labelgb <- factor(c(rep("ALL",27),rep("AML",11)))
names(golub) <- labelgb
golub2<-log2(golub)
plot(golub2)
biocLite('limma')
library('limma')
fit3 = lmFit(golub2)
e3 = eBayes(fit3)
tab3<-topTable(e3, number=nrow(golub))
head(tab3)
fit3@.Data[[9]]
fit4 = lmFit(golub2, design)
```


```{r shuffle}
ecoli <- read.delim("Y:/felles/PROTEOMICS and XRAY/Ani/Christina/Strainssamples.txt", header=TRUE, na.strings='', sep="\t")
for  (cnt in 1:6){
  cnt
}
  ecoli_shuffle <- ecoli_shuffle[sample(nrow(ecoli)),]
write.csv(file="Y:/felles/PROTEOMICS and XRAY/Ani/Christina/Strainssamples_shuffled_rn.csv",x=ecoli_shuffle,row.names = F)
```


```{r prot}
install.packages("PTXQC")
library("PTXQC")
help(package = "PTXQC")
vignette("PTXQC-InputData", package = "PTXQC")
txt_folder = "F:/promec/Elite/LARS/2016/desember/Lymphoma/Results"
r = createReport(txt_folder)
cat(paste0("\nReport generated as '", r$report_file, "'\n\n"))

```


```{r Cluster}
colscl <- colorRampPalette(c("darkblue", "chocolate"), space = "rgb")(100)
heatmap(na.omit(as.matrix(data)),col=colscl)
```




```{r SVA}
#https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btu375
#source("http://www.bioconductor.org/biocLite.R")
#biocLite(c("Biobase","sva","bladderbatch","snpStats"))
  library(sva)
  library(devtools)
  library(Biobase)
  library(bladderbatch)
  library(snpStats)

data(bladderdata)
pheno = pData(bladderEset)
edata = exprs(bladderEset)

mod = model.matrix(~as.factor(cancer) + as.factor(batch),data=pheno)
fit = lm.fit(mod,t(edata))
hist(fit$coefficients[2,],col=2,breaks=100)
table(pheno$cancer,pheno$batch)

batch = pheno$batch
modcombat = model.matrix(~1, data=pheno)
modcancer = model.matrix(~cancer, data=pheno)
combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)

combat_fit = lm.fit(modcancer,t(combat_edata))
hist(combat_fit$coefficients[2,],col=2,breaks=100)

plot(fit$coefficients[2,],combat_fit$coefficients[2,],col=2,
      xlab="Linear Model",ylab="Combat",xlim=c(-5,5),ylim=c(-5,5))
abline(c(0,1),col=1,lwd=3)

mod = model.matrix(~cancer,data=pheno)
mod0 = model.matrix(~1, data=pheno)
sva1 = sva(edata,mod,mod0,n.sv=2)

summary(lm(sva1$sv ~ pheno$batch))

boxplot(sva1$sv[,2] ~ pheno$batch)
points(sva1$sv[,2] ~ jitter(as.numeric(pheno$batch)),col=as.numeric(pheno$batch))

modsv = cbind(mod,sva1$sv)
fitsv = lm.fit(modsv,t(edata))

par(mfrow=c(1,2))
plot(fitsv$coefficients[2,],combat_fit$coefficients[2,],col=2,
      xlab="SVA",ylab="Combat",xlim=c(-5,5),ylim=c(-5,5))
abline(c(0,1),col=1,lwd=3)
plot(fitsv$coefficients[2,], fit$coefficients[2,],col=2,
      xlab="SVA",ylab="linear model",xlim=c(-5,5),ylim=c(-5,5))
abline(c(0,1),col=1,lwd=3)

data(for.exercise)
controls <- rownames(subject.support)[subject.support$cc==0]
use <- seq(1, ncol(snps.10), 10)
ctl.10 <- snps.10[controls,use]

xxmat <- xxt(ctl.10, correct.for.missing=FALSE)
evv <- eigen(xxmat, symmetric=TRUE)
pcs <- evv$vectors[,1:5]

pop <- subject.support[controls,"stratum"]
plot(pcs[,1],pcs[,2],col=as.numeric(pop),
      xlab="PC1",ylab="PC2")
legend(0,0.15,legend=levels(pop),pch=19,col=1:2)
```

---
title: "batch-effect-analysis"
author: "Ani"
date: "Monday, Feb 15, 2016"
output: html_document
---