#####********************************** MOA sRNAomics- Script to determine statistically significant differences in sRNA levels**********************************######
# Last modified: 2019-06-20

#####************************************************* Workspace Preparation ************************************************#####

# Working directory:
setwd("..")

# Required libraries


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("limma")

library(limma)


if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("edgeR")

library (edgeR)

#####**************************************************Data analysis********************************************************#####

### Step 1:  Preparation and Import of the sRNA read count data

# miRNA and phasiRNA read counts obtained by sRNA-Seq are given in MOA2010-05_3_sRNAomics_sRNA-Seq_miRNAs_counts.xlsx and MOA2010-05_3_sRNAomics_sRNA-Seq_phasiRNAs_counts.xlsx, respectively.
# Preparation of count table (i.e. count_table.txt) by combining miRNA and phasiRNA counts from these two files. 
# Count table should contain sRNAs in rows (first column name is sRNAs), and sRNA read counts in columns (column names are the same as sample IDs). 
# Sequence of the samples as follows: MOAs_1269,MOAs_1271,MOAs_1273,MOAs_1562,MOAs_1564,MOAs_1566,MOAs_0666,MOAs_0668,MOAs_0670,MOAs_0961,MOAs_0963,MOAs_0965

# Reading in count table into object x :

x <- read.table("count_table.txt",header=TRUE,row.names="sRNAs",sep="\t", stringsAsFactors=FALSE)

#####********************************************** sRNA differential expression analysis ***********************************************#####

### Step 2.1. Defining groups
group <- factor(c(1,1,1,2,2,2,3,3,3,4,4,4))

### Step 2.2. Data conversion into a DGEList-object using the DGEList function and specification of the library sizes of the samples.
y <- DGEList(counts=x,group=group,
lib.size=c(2866133,3214236,2592971,4016294,3094686,3936757,1191167,2818468,2005330,2252675,3727295,1961316))

### Step 2.3. Conversion of raw counts to CPM
cpm <- cpm(y)

### Step 2.4. Removing sRNAs that are lowly expressed – only sRNAs having CPM value > 1 in at least three samples are kept.
keep.exprs <- rowSums(cpm>1)>=3

### Step 2.5. Data normalization
y1 <- y[keep.exprs, , keep.lib.sizes=TRUE]
y1 <- calcNormFactors(y1)

### Step 2.5.Creating design matrix and contrasts
design = cbind(Desiree_mock = c(1,1,1,0,0,0,0,0,0,0,0,0), Desiree_PVY = c(0,0,0,1,1,1,0,0,0,0,0,0), NahG_mock =
c(0,0,0,0,0,0,1,1,1,0,0,0), NahG_PVY = c(0,0,0,0,0,0,0,0,0,1,1,1))
contrastMatrix = makeContrasts("Desiree_PVY-Desiree_mock","NahG_PVY-NahG_mock","NahG_mock-Desiree_mock","NahG_PVY-Desiree_
PVY", levels=design)

### Step 2.5. Data transformation by voom
v <- voom(y1,design,plot=TRUE)

### Step 2.6. Fitting linear models
fit <- lmFit(v,design)
fit2 = contrasts.fit(fit, contrastMatrix)
fit2 <- eBayes(fit2)

### Step 2.7. Outputting the results
Desiree_PVY_mock <- topTable(fit2, coef=1, number=1000000)
write.table(Desiree_PVY_mock,file= "Desiree_PVY_mock.txt", sep="\t")
NahG_PVY_mock <- topTable(fit2, coef=2, number=1000000)
write.table(NahG_PVY_mock,file= "NahG_PVY_mock.txt", sep="\t")
NahG_mock_Desiree_mock <- topTable(fit2, coef=3, number=1000000)
write.table(NahG_mock_Desiree_mock,file= "NahG_mock_Desiree_mock.txt", sep="\t")
NahG_PVY_Desiree_PVY <- topTable(fit2, coef=4, number=1000000)
write.table(NahG_PVY_Desiree_PVY,file= "NahG_PVY_Desiree_PVY.txt", sep="\t")