Scripts

###### Systems biology analysis to understand regulatory miRNA networks in lung cancer

#Kunz M, Pittroff A, Dandekar T. Systems biology analysis to understand 
# regulatory miRNA networks in lung cancer. Methods Mol Biol. 2018, vol. 
# 1819, https://doi.org/10.1007/978-1-4939-8618-7_11
#3.1 Software and package installation (code requires R version 3.2.2)

source("http://bioconductor.org/biocLite.R")
biocLite()
biocLite("GEOquery")
biocLite("limma")
biocLite("gplots")
#3.2 Load the R libraries
library("GEOquery")
library("limma")
library("gplots")
#3.3 Downloading the miRNA expression dataset from GEO
GSE63805 <- getGEO("GSE63805")	# GEO ID of interest can be used
eset1 = GSE63805[[1]]
eset1b=eset1
#3.4 Data preparation
colnames(pData(eset1b))[8] = "Tumor"
levels(eset1b$Tumor) = c("Control","Lung Cancer")
dim(eset1b)
length(eset1b$Tumor[eset1b$Tumor == "Control"])
length(eset1b$Tumor[eset1b$Tumor == "Lung Cancer"])
#3.5 Check Normalization
boxplot(exprs(eset1b))
ctrl<-which(eset1b$Tumor=="Control")[1]
test<-which(eset1b$Tumor=="Lung Cancer")[1]
exp<-exprs(eset1b)[,c(ctrl,test)]
plot(exp)
exprs(eset1b)=log2(exprs(eset1b))
#Visualize the log2 transformed expression values: repeat steps 3.5 1 to 5 to generate the plots
#3.6 Selecting differentially expressed genes
labels=pData(eset1b)[ ,"Tumor"]
design=model.matrix(~labels)
fit=eBayes(lmFit(eset1b, design))
t<-topTable(fit,coef=2,number=Inf,p.value=0.05,lfc=2)	# p.value sets a cutoff for adjusted p-values, lfc sets a cutoff for minimum logFC
l.deg=list()
l.deg<-intersect(rownames(exprs(eset1b)),rownames(t))
write.table(l.deg, file="topDEGs.txt", quote = FALSE, sep = "\t", row.names=FALSE, col.names=FALSE)
#3.7 Data Visualization of differentially expressed genes (see Note 4)
color.map<-function(Tumor){if(Tumor=="Control") "#0000FF" else "#FF0000"}    # adding colorcode for sample-identification (control=blue, tumor=red)
patientcolors<-unlist(lapply(eset1b$Tumor, color.map))
heatmap.2(exprs(eset1b[l.deg,]),ColSideColors=patientcolors, key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.8)
############################################################