My_Code_Transcriptomics

library(Rsubread)
library(edgeR)
library(limma)
library(org.Mm.eg.db)
library(annotate)

# Building the index in order to speed up the queries
buildindex(basename="Indexed_SmallGenomeMouse",reference="SmallGenomeMouse.fa.gz")

### This line of code indicates to the Rsubread aligner to only consider the fastq.gz files.  
all_12_fastq <- list.files(path = "C:/Users/blado/Desktop/THE ANALYSIS/ReferenceGenome", 
                          pattern = ".fastq.gz$", full.names = TRUE)
### All the fastq.gz files containing the reads are now in this vector. 
all_12_fastq
#### All the  fastq.gz will be now aligned to the index that we have built
align(index = "Indexed_SmallGenomeMouse", readfile1 = all_12_fastq)


## This line of code will cluster all the generated Bam files and put them into this vector. 
all_12_bam <- list.files(path = "C:/Users/blado/Desktop/THE ANALYSIS/ReferenceGenome",
                        pattern = ".BAM$", full.names = TRUE)
all_12_bam
## Recapitulation of the amount of reads that mapped to SmallGenomeMouse
Prop <- propmapped(files=all_12_bam)
Prop


## Let's use the featureCounts function to count the mapped reads 
the_count <- featureCounts(all_12_bam, annot.inbuilt="mm10")
# Checking the slots stored in the_count
names(the_count)
## Let's peek in the stats 
the_count$stat
## Let's peek in the dimensions to observe the number of genes
dim(the_count$counts)
##  These are the first 6 lines
head(the_count$counts)
head(the_count$annotation)


#### Let's write counts to table 
the_count_counts <- the_count$counts
### checking column names
colnames(the_count_counts)
#### changing columns names
colnames(the_count_counts) <- c("DG-BV", "DH-BV", "DI-BP", "DJ-BP", "DK-BL",
                                "DL-BL", "LA-LV", "LB-LV", "LC-LP", "LD-LP",
                                "LE-LL", "LF-LL")
### Converting the_count_counts into a data frame
the_count_counts <- as.data.frame(the_count_counts)
### Let's sort the data frame 
Sorted_counts <- the_count_counts[order(-the_count_counts$`DG-BV`),]
### Converting my row names to a proper column and renaming it. 
Sorted_counts <- cbind(GeneID = rownames(Sorted_counts), Sorted_counts)
### Removing the original row names
rownames(Sorted_counts) <- NULL

### Let's convert the list of GeneID into Gene Name with the getSYMBOL function
Gene_Name <- getSYMBOL(Sorted_counts$GeneID, data='org.Mm.eg.db')
### Merging the obtained gene name with the data frame "Sorted_counts"
Sorted_counts1 <- cbind(Gene_Name, Sorted_counts)
### Removing the again the original row names
rownames(Sorted_counts1) <- NULL

### Let's output the count table in csv format
write.csv(x=Sorted_counts1, file="COUNT_TABLE.csv",row.names=TRUE)
#Let's load data 
CountTable <- read.csv(file = 'COUNT_TABLE.csv', header = TRUE)

# Let's Read the sample information
### We first converted SampleTable.csv into SampleTable.txt
### Then we have renamed it into SampleInfo.txt
### And finally, we have rearranged the information inside the file 
sampleinfo <- read.delim("SampleInfo.txt", header=TRUE, sep=",")
View(sampleinfo)


### Let's read the full-count-table in csv format into R
fullcounttable <- read.csv("full-count-table.csv", sep = ",", header = TRUE)
#### changing columns names
colnames(fullcounttable ) <- c("EntrezGeneID","Length","DG-BV", "DH-BV", "DI-BP", "DJ-BP", "DK-BL",
                                "DL-BL", "LA-LV", "LB-LV", "LC-LP", "LD-LP",
                                "LE-LL", "LF-LL")
### convert the factors to strings/characters 
fullcounttable$EntrezGeneID <- as.character(fullcounttable$EntrezGeneID)
### Let's convert the EntrezGeneID into Gene Symbol with the getSYMBOL function
GeneSymbol <- getSYMBOL(fullcounttable$EntrezGeneID, data='org.Mm.eg.db')
### Merging Gene_Symbol with fullcounttable
fullcounttable <- cbind(GeneSymbol, fullcounttable)
head(fullcounttable)
###  Let's remove first three columns
Count_Data <- fullcounttable[,-(1:3)]
# Let's store EntrezGeneID as rownames
rownames(Count_Data) <- fullcounttable[,2]
## Viewing the rownames
rownames(Count_Data)


#### FILTERING LOWLY EXPRESSED GENES

# Getting the CPMs
My_CPM <- cpm(Count_Data)
# Let's peek in the output
head(My_CPM)
# Selecting the values in My_CPM that are greater than 0.5
Thresh_val <- My_CPM > 0.5
# Yielding the logical matrix with TRUEs and FALSEs
head(Thresh_val)


# Recapitulation of the number of TRUEs in each row
# 11433 genes having TRUEs in all 12 samples.
table(rowSums(Thresh_val))
# Let's keep genes that have at least 2 TRUES in each row 
keeping <- rowSums(Thresh_val) >= 2

# A way of retaining the more highly expressed genes in the rows of Count_Data
Counts_keep <- Count_Data[keeping,]
summary(keeping)
## Mode   FALSE    TRUE
## logical   11375   15804 
dim(Counts_keep)
# [1] 15804    12


# Looking at the first sample
plot(My_CPM[,1],Count_Data[,1], main=colnames(My_CPM)[1])


# Limiting the x and y-axis to see what is happening at the smaller counts
plot(My_CPM[,1],Count_Data[,1],ylim=c(0,50),xlim=c(0,3), main=colnames(My_CPM)[1])
#Placing Vertical line at 0.5 CPM
abline(v=0.5)


### from counts to DGEList object
Dge_Obj <- DGEList(Counts_keep)
### Let's view the created Dge_Obj
Dge_Obj
# Slots stored in Dge_Obj
names(Dge_Obj)
# In the samples slot we find library size information
Dge_Obj$samples
## Getting the number of reads for each sample
Dge_Obj$samples$lib.size


## Barplot Generation of the library sizes
barplot(Dge_Obj$samples$lib.size, names=colnames(Dge_Obj), las=2)
# Let's entitle the plot
title("Barplot of library sizes")

### Let's normalise the DGEList object
Dge_Obj_Nor <- calcNormFactors(Dge_Obj)
### Let's peek in the new normalisation factors
Dge_Obj_Nor$samples


### Treatment factors are combined into a unique grouping factor
Group <- paste(sampleinfo$CellType, sampleinfo$Status, sep=".")
Group <- factor(Group)

### Categorical tabulation of the factors 
table(Group)

### MDS plot decoration step
pch <- c(0,1,2,15,16,17)
colors <- rep(c("darkgreen", "red", "blue"), 2)
## Creating another plotting window
dev.new(width=5, height=4)
## Clustering in two dimensions the samples with MDS plot
plotMDS(Dge_Obj, col=colors[Group], pch=pch[Group])
## Adding a legend to the MDS plot
legend("topleft", legend=levels(Group), pch=pch, col=colors, ncol=2)


### Step1: Pre-Normalization
# log2 counts per million before normalization
logcount_UN <- cpm(Dge_Obj,log=TRUE)
# Verifying the distributions of samples with boxplots
boxplot(logcount_UN, xlab="", ylab="Log2 counts per million",las=2)
# Adding a red horizontal line which is the median logCPM
abline(h=median(logcount_UN),col="red")
## Let's title the plot
title("Boxplots of logCPMs (unnormalised)")

### Step2:  Post-Normalization
# log2 counts per million
logcount_Nor <- cpm(Dge_Obj_Nor,log=TRUE)
# Verifying the distributions of samples with boxplots
boxplot(logcount_Nor, xlab="", ylab="Log2 counts per million",las=2)
#  Adding a red horizontal line which is the median logCPM
abline(h=median(logcount_Nor),col="red")
### Let's title the plot
title("Boxplots of logCPMs (Normalised)")


### The mean-difference plots showing the composition bias problem
### MD plot Pre-TMM normalization
plotMD(Dge_Obj,column = 7)
abline(h=0,col="green", lwd = 2, lty = 2)

### MD plot Post-TMM normalization
plotMD(Dge_Obj_Nor,column = 7)
abline(h=0,col="green", lwd = 2, lty = 2)






############ Differential expression analysis of "basal pregnant"
############  vs "basallactate"

### design matrix construction
Design <- model.matrix(~0+Group)
### Naming the Design columns with the levels names
colnames(Design) <- levels(Group)
### Printing the design matrix
Design



### Dispersion estimation 
Dge_Obj_Dis <- estimateDisp(Dge_Obj_Nor, Design, robust=TRUE)
### Viewing the slots stored in Dge_Obj_Dis
names(Dge_Obj_Dis)
### Dispersion estimation vizualization with plotBCV
plotBCV(Dge_Obj_Dis)

### Let's use the glmQLFit function to estimate the QL dispersions
Fit <- glmQLFit(Dge_Obj_Dis, Design, robust=TRUE)
### Let us Plot the QL Dispersions
plotQLDisp(Fit)





### Construction of contrast
BP_BL <- makeContrasts(basal.pregnant-basal.lactate, levels = Design) 
### Performing a QL F-test
Res <- glmQLFTest(Fit, contrast = BP_BL)
### Printing the dimension of "Res"
dim(Res) ### [1] 15804     4
### Using the topTags function to extract the top differentially expressed genes
Top_Tags <- topTags(Res,n=15804)
### Printing the five first rows of the object of class "TopTags"  
head(Top_Tags)
# Let's take out the data.frame incorporating DE results for the top genes
Top_Tags_Table <- Top_Tags$table
### row names represents Gene ID
Gene_ID = rownames(Top_Tags_Table)
### Creating a new Gene_ID column
Top_Tags_Table <- cbind(Gene_ID, Top_Tags_Table)
### Converting the list of Gene ID into Gene Symbol with the getSYMBOL function
Gene_Symbol <- getSYMBOL(Top_Tags_Table$Gene_ID, data='org.Mm.eg.db')
### Merging Gene_Symbol with Top_Tags_Table
Top_Tags_Table <- cbind(Gene_Symbol, Top_Tags_Table)




### Determination of the total amount of DE genes 
is.de <- decideTestsDGE(Res, p.value=0.05)
### Recapitulation
summary(is.de)
##       -1*basal.lactate 1*basal.pregnant
## Down                                2505
## NotSig                             10529
## Up                                  2770
### Picturing the extent of the DE changes
plotMD(Res, status=is.de, values=c(1,-1), col=c("red","blue"),
       legend="topright")





#### DE above a fold-change threshold
### Investigation of DE fold changes higher than 1.5
TR <- glmTreat(Fit, contrast=BP_BL, lfc=log2(1.5))
### Printing the dimension of "TR"
dim(TR) ### [1] 15804     4
### topTags function for extraction of DE genes
Top_Tags1 <- topTags(TR,n=15804)
### Printing the five first rows of the object of class "TopTags"
head(Top_Tags1)
### Let's bring out the data.frame incorporating DE results for the top genes
Top_Tags_Table1 <- Top_Tags1$table
### row names represents Gene ID
Gene_ID = rownames(Top_Tags_Table1)
### Creating a new Gene_ID column
Top_Tags_Table1 <- cbind(Gene_ID, Top_Tags_Table1)
### Converting the list of Gene ID into Gene Symbol with the getSYMBOL function
Gene_Symbol <- getSYMBOL(Top_Tags_Table1$Gene_ID, data='org.Mm.eg.db')
### Merging Gene_Symbol with Top_Tags_Table1
Top_Tags_Table1 <- cbind(Gene_Symbol, Top_Tags_Table1)




### Finding the total amount of DE genes
is.de1 <- decideTestsDGE(TR, p.value=0.05)
### Recapitulation 
summary(is.de1)
##       -1*basal.lactate 1*basal.pregnant
## Down                                 366
## NotSig                             14761
## Up                                   677
### Results visualization 
plotMD(TR, status=is.de1, values=c(1,-1), col=c("red","blue"),
       legend="topright")


############ Differential expression analysis of Basal vs Luminal



### Creating a design matrix

## Extracting the treatment factor CellType
Group1 <- paste(sampleinfo$CellType)
## Creating the factor
Group1 <- factor(Group1)
## Performing a categorical tabulation of the factors
table(Group1) ## basal luminal
              ##     6       6

## Let us construct the design matrix without intercept
Design1 <- model.matrix(~0+Group1)
## Defining new col names
colnames(Design1) <- levels(Group1)
### Printing "Design1"
Design1



### Dispersion estimation
Dge_Obj_Dis1 <- estimateDisp(Dge_Obj_Nor, Design1, robust=FALSE)
### slots stored in Dge_Obj_Dis1 
names(Dge_Obj_Dis1)
### Viewing the dispersion with the plotBCV
plotBCV(Dge_Obj_Dis1)



### QL dispersion estimation
Fit1 <- glmQLFit(Dge_Obj_Dis1, Design1, robust=TRUE)
### Plotting the QL Dispersion
plotQLDisp(Fit1)



### Building the contrast
B_L <- makeContrasts(basal-luminal, levels = Design1) 
### QL F-test
Res1 <- glmQLFTest(Fit1, contrast = B_L)
### Dimension of "Res1"
dim(Res1) ### [1] 15804     4
### Extraction of DE genes
Top_Tags2 <- topTags(Res1,n=15804)
### Let's peek in "Top_Tags2"  
head(Top_Tags2)
# Data frame with DE genes
Top_Tags_Table2 <- Top_Tags2$table

#### Adding annotation
### row names = Gene ID
Gene_ID = rownames(Top_Tags_Table2)
### New Gene_ID column
Top_Tags_Table2 <- cbind(Gene_ID, Top_Tags_Table2)
### From Gene ID to Gene Symbol
Gene_Symbol <- getSYMBOL(Top_Tags_Table2$Gene_ID, data='org.Mm.eg.db')
### Merging Gene_Symbol with Top_Tags_Table
Top_Tags_Table2 <- cbind(Gene_Symbol, Top_Tags_Table2)




### The total amount of DE genes 
is.de2 <- decideTestsDGE(Res1, p.value=0.05)
### Summarization
summary(is.de2)
##       1*basal -1*luminal
##Down                 4783
##NotSig               6233
##Up                   4788 
### Magnitude of the DE changes
plotMD(Res1, status=is.de2, values=c(1,-1), col=c("red","blue"),
       legend="topright")




### DE above a fold change threshold (Basal vs Luminal)

### DE fold changes higher than 3
TR1 <- glmTreat(Fit1, contrast=B_L, lfc=log2(3))
### dim of "TR1"
dim(TR1) ### [1] 15804     4
### extraction of DE genes
Top_Tags3 <- topTags(TR1,n=15804)
### Displaying first rows
head(Top_Tags3)
### Table with DE genes
Top_Tags_Table3 <- Top_Tags3$table
### Let us add Gene Symbol
## Row names are Gene ID
Gene_ID = rownames(Top_Tags_Table3)
###  New Gene_ID column
Top_Tags_Table3 <- cbind(Gene_ID, Top_Tags_Table3)
### Converting 
Gene_Symbol <- getSYMBOL(Top_Tags_Table3$Gene_ID, data='org.Mm.eg.db')
### Merging 
Top_Tags_Table3 <- cbind(Gene_Symbol, Top_Tags_Table3)




### Global amount of DE genes
is.de3 <- decideTestsDGE(TR1, p.value=0.05)
### Recap
summary(is.de3)
##       1*basal -1*luminal
## Down                 972
## NotSig             13104
## Up                  1728
### visualization 
plotMD(TR1, status=is.de3, values=c(1,-1), col=c("red","blue"),
       legend="topright")