FM Figures about RNAseq Analaysis - Rachelle
Gene expression along the genome
#Set the "RNA_seq_data" directory as working directory.
setwd("~/Library/CloudStorage/OneDrive-UniversitéParis-Saclay/Fundamental Microbiology (2026-2027)/Functional Genomics/RNA_seq-20261022") #Replace by your path
RNA<-read.table("tables/MP48vsMP24.complete.txt", header=T) #We read the table
summary(RNA) #We have an quick overview of the content of this table## Id MP24a MP24b MP48a
## Length :7534 Min. : 0 Min. : 0 Min. : 0.0
## N.unique :7534 1st Qu.: 44 1st Qu.: 37 1st Qu.: 44.2
## N.blank : 0 Median : 290 Median : 242 Median : 211.0
## Min.nchar: 18 Mean : 3520 Mean : 3101 Mean : 2536.7
## Max.nchar: 21 3rd Qu.: 1534 3rd Qu.: 1343 3rd Qu.: 984.0
## Max. :3026484 Max. :3032883 Max. :3345901.0
##
## MP48b Ycongo24a Ycongo24b Y24a
## Min. : 0 Min. : 0.0 Min. : 0.0 Min. : 0
## 1st Qu.: 48 1st Qu.: 24.0 1st Qu.: 18.0 1st Qu.: 26
## Median : 235 Median : 159.5 Median : 123.0 Median : 187
## Mean : 2573 Mean : 2710.6 Mean : 1683.2 Mean : 2524
## 3rd Qu.: 1031 3rd Qu.: 851.8 3rd Qu.: 648.8 3rd Qu.: 1097
## Max. :5619211 Max. :6811646.0 Max. :3655746.0 Max. :3106591
##
## Y24b Y48a Y48b norm.MP24a
## Min. : 0.0 Min. : 0.0 Min. : 0 Min. : 0
## 1st Qu.: 36.0 1st Qu.: 37.0 1st Qu.: 54 1st Qu.: 29
## Median : 239.5 Median : 175.0 Median : 265 Median : 193
## Mean : 2474.0 Mean : 1608.3 Mean : 2233 Mean : 2340
## 3rd Qu.: 1262.0 3rd Qu.: 747.8 3rd Qu.: 1106 3rd Qu.: 1020
## Max. :2470345.0 Max. :2315818.0 Max. :3463663 Max. :2012191
##
## norm.MP24b norm.MP48a norm.MP48b norm.Ycongo24a
## Min. : 0 Min. : 0.0 Min. : 0 Min. : 0.0
## 1st Qu.: 29 1st Qu.: 47.2 1st Qu.: 50 1st Qu.: 29.0
## Median : 191 Median : 226.0 Median : 245 Median : 195.5
## Mean : 2447 Mean : 2721.9 Mean : 2681 Mean : 3321.9
## 3rd Qu.: 1060 3rd Qu.: 1056.0 3rd Qu.: 1074 3rd Qu.: 1043.8
## Max. :2393442 Max. :3590152.0 Max. :5855550 Max. :8347824.0
##
## norm.Ycongo24b norm.Y24a norm.Y24b norm.Y48a
## Min. : 0 Min. : 0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 28 1st Qu.: 27 1st Qu.: 29.0 1st Qu.: 49.0
## Median : 194 Median : 194 Median : 195.5 Median : 232.0
## Mean : 2656 Mean : 2621 Mean : 2019.9 Mean : 2128.2
## 3rd Qu.: 1024 3rd Qu.: 1139 3rd Qu.: 1030.0 3rd Qu.: 989.5
## Max. :5767543 Max. :3225987 Max. :2016963.0 Max. :3064405.0
##
## norm.Y48b baseMean MP24 MP48
## Min. : 0.0 Min. : 0.0 Min. : 0 Min. : 0
## 1st Qu.: 49.0 1st Qu.: 60.1 1st Qu.: 30 1st Qu.: 50
## Median : 239.0 Median : 313.5 Median : 196 Median : 240
## Mean : 2017.6 Mean : 2495.5 Mean : 2394 Mean : 2702
## 3rd Qu.: 998.8 3rd Qu.: 1285.0 3rd Qu.: 1060 3rd Qu.: 1076
## Max. :3129447.0 Max. :3940350.5 Max. :2202816 Max. :4722851
##
## Ycongo24 Y24 Y48 FoldChange
## Min. : 0.0 Min. : 0 Min. : 0.0 Min. : 0.005
## 1st Qu.: 30.0 1st Qu.: 29 1st Qu.: 50.0 1st Qu.: 0.477
## Median : 197.5 Median : 198 Median : 236.0 Median : 1.013
## Mean : 2988.7 Mean : 2321 Mean : 2072.9 Mean : 12.858
## 3rd Qu.: 1048.0 3rd Qu.: 1094 3rd Qu.: 993.5 3rd Qu.: 2.654
## Max. :7057684.0 Max. :2621475 Max. :3096926.0 Max. :3130.961
## NAs :213
## log2FoldChange stat pvalue padj
## Min. :-7.6040 Min. :-22.5880 Min. :0.000000 Min. :0.000000
## 1st Qu.:-1.0690 1st Qu.: -3.1570 1st Qu.:0.000000 1st Qu.:0.000000
## Median : 0.0190 Median : 0.0420 Median :0.001212 Median :0.002424
## Mean : 0.3439 Mean : 0.6696 Mean :0.143914 Mean :0.162861
## 3rd Qu.: 1.4080 3rd Qu.: 3.3740 3rd Qu.:0.158841 3rd Qu.:0.211778
## Max. :11.6120 Max. : 40.6180 Max. :1.000000 Max. :1.000000
## NAs :213 NAs :213 NAs :213 NAs :213
## dispGeneEst dispFit dispMAP dispersion
## Min. : 0.00000 Min. : 0.0474 Min. : 0.0147 Min. : 0.0147
## 1st Qu.: 0.00000 1st Qu.: 0.0492 1st Qu.: 0.0335 1st Qu.: 0.0335
## Median : 0.02100 Median : 0.0544 Median : 0.0498 Median : 0.0498
## Mean : 0.05693 Mean : 0.1667 Mean : 0.1394 Mean : 0.1443
## 3rd Qu.: 0.05460 3rd Qu.: 0.0806 3rd Qu.: 0.0858 3rd Qu.: 0.0858
## Max. :10.00000 Max. :35.7365 Max. :10.0000 Max. :10.0000
## NAs :213 NAs :213 NAs :213 NAs :213
## betaConv maxCooks
## Mode:logical Mode:logical
## TRUE:7321 NAs :7534
## NAs :213
##
##
##
## #We need to annotate the file to indicate the replicon and start-stop positions of each
#Set the "Counts" directory as working directory, because we need also this information that comes from another table
ANNOT<-read.table("Counts/WTMP5_24_1_counts_S.txt", header=T)
#We conserve only the 6 first columns
ANNOT<-ANNOT[,-c(7)]
#We rename the column (so that the "Id" will be present in the annotation file as well as our RNA file, that's the common column)
colnames(ANNOT)<-c("Id", "Replicon", "START","End", "Strand", "Length")
#We fuse the tables:
RNA<-merge(RNA, ANNOT, by="Id", all.x=T)
#If you want you can save this table with this command line
write.csv2(RNA, "RNA.csv", row.names = F)
# We will analyzed only the chromosomal genes (without genes from the plasmid). To do this, we can use the "subset" function like this:
RNA_CHR<-subset(RNA, RNA$Replicon=="NZ_CP012382.1")
#Order RNA_CHR by the start position, first one the lowest, latest the
RNA_CHR$START<-as.numeric(RNA_CHR$START)
RNA_CHR<-RNA_CHR[order(RNA_CHR$START),]#We order the rows depending on the start value
#If you want you can save this table with this command line
write.csv2(RNA_CHR, "RNA_CHR.csv", row.names = F)
# Calculation of the mean sense transcription in a window of 81 genes:
RNA_CHR$RNAseq_24h_WINDOW<-"NA" # This way, we create a new variable in which we will put the information concerning the mean persistence in a window centered at each gene position
RNA_CHR$RNAseq_48h_WINDOW<-"NA"
a<-20 #size of the window, before and after each gene position (changed to 20)
b<-a+1 #the calculation will start at this position in the genome of interest (this information will change automatically of we change the size of the window)
l<-length(RNA_CHR$Id)-a #last position analyzed in the genome (=> number of genes minus the half size of the windows, given that the windows are centered at each gene position)
for (i in b:l){
j = i - a
k = i + a
RNA_CHR$RNAseq_24h_WINDOW[i]<-mean(RNA_CHR$MP24[j:k], na.rm=TRUE)
RNA_CHR$RNAseq_48h_WINDOW[i]<-mean(RNA_CHR$MP48[j:k], na.rm=TRUE)
}
#I always convert the outcome of a loop as numeric values (since sometimes you can have troubles with values that could not be calculated, etc...)
RNA_CHR$RNAseq_24h_WINDOW<-as.numeric(RNA_CHR$RNAseq_24h_WINDOW)
RNA_CHR$RNAseq_48h_WINDOW<-as.numeric(RNA_CHR$RNAseq_48h_WINDOW)
# Plot #1 - sense transcription along the chromosome in exponential phase (type = "l" is indicated to represent a line instead of dots):
plot(RNA_CHR$START, RNA_CHR$RNAseq_24h_WINDOW, type = "l", xlab ="Position along the chromosome (bp)", ylab = "Mean transcription (in the window)")
# Plot #2 - sense transcription along the chromosome in stationary phase:
plot(RNA_CHR$START, RNA_CHR$RNAseq_48h_WINDOW, type = "l", xlab ="Position along the chromosome (bp)", ylab = "Mean transcription (in the window)")
# Plot #2 - sense transcription along the chromosome in stationary phase:
plot(RNA_CHR$START, RNA_CHR$RNAseq_48h_WINDOW, type = "l", xlab ="Position along the chromosome (bp)", ylab = "Mean transcription (in the window)")
# Plot #3 - overlay:
plot(RNA_CHR$START, RNA_CHR$RNAseq_24h_WINDOW, type = "l", xlab ="Position along the chromosome (bp)", ylab = "Mean transcription (in the window)")
lines(RNA_CHR$START, RNA_CHR$RNAseq_48h_WINDOW, col = "pink")
# We can also add information concerning the position of the origin of replication (in red), of the ribosomal RNA encoding operons (in blue) and of the TIR (terminal inverted repeats, in grey), like this:
abline(v=c(196959, 4023406, 1473191, 1920301, 3247066, 3784403, 4485769, 5802474), col=c("grey","red", "blue","blue","blue","blue","blue","blue"), lty=2, lwd=1)
# Let's start by calculating the log2 values and redo like above. I have to be carefull with the 0... (that will induce a bug with the log2 function ;-))
RNA_CHR$RNAseq_24h_log2<-"NA"
RNA_CHR$RNAseq_48h_log2<-"NA"
#To avoid the "zero problem", we add "1" to the number of reads in each condition
RNA_CHR$RNAseq_24h_log2<-log2(RNA_CHR$MP24+1)
RNA_CHR$RNAseq_48h_log2<-log2(RNA_CHR$MP48+1)
RNA_CHR$RNAseq_24h_WINDOW_log2<-"NA"
RNA_CHR$RNAseq_48h_WINDOW_log2<-"NA"
a<-20
b<-a+1
l<-length(RNA_CHR$Id)-a
for (i in b:l){
j = i - a
k = i + a
RNA_CHR$RNAseq_24h_WINDOW_log2[i]<-mean(RNA_CHR$RNAseq_24h_log2[j:k], na.rm=TRUE)
RNA_CHR$RNAseq_48h_WINDOW_log2[i]<-mean(RNA_CHR$RNAseq_48h_log2[j:k], na.rm=TRUE)
}
RNA_CHR$RNAseq_24h_WINDOW_log2<-as.numeric(RNA_CHR$RNAseq_24h_WINDOW_log2)
RNA_CHR$RNAseq_48h_WINDOW_log2<-as.numeric(RNA_CHR$RNAseq_48h_WINDOW_log2)
plot(RNA_CHR$START, RNA_CHR$RNAseq_24h_WINDOW_log2, type = "l", xlab ="Position along the chromosome (bp)", ylab = "Mean transcription (in the window, log2)")
lines(RNA_CHR$START, RNA_CHR$RNAseq_48h_WINDOW_log2, col = "purple")
abline(v=c(196959, 4023406, 1473191, 1920301, 3247066, 3784403, 4485769, 5802474), col=c("grey","red", "blue","blue","blue","blue","blue","blue"), lty=2, lwd=1)
Analyze the figures. What biological phenomena do you observe? What is
the interest of representing logged values? R/ we can observe in the
last figure that Streptomyces has consistent high expression of genes
located in the central core of the genome, where the origin of
replication (Ori) is located. This shows active replication and
transcription. We also observe high expression of ribosomal genes. And
then at the peripheral regions the gene expression decrease due to the
presence of accesory gene clusters. The logged values help to visualize
the data in a clearer way.
Visualization of the data using heatmap representation
We will work with the same RNA_CHR table that we have previously generated. First, we will normalize the mean DESeq2 normalized counts on gene size. We will obtain normRPK (normalized reads per kb). Then we will perform a heatmap representation of the data using the function heatmap(). If you indicate, ‘?heatmap’ in the terminal you can have detailed information about this function. If you indicate just ‘heatmap’, you will have its code.
#Read the RNA_CHR tabvle that you have previously generated
RNA_CHR<-read.csv2("RNA_CHR.csv")
# First create new columns containing the normRPK value
RNA_CHR$RPK_MP24<-RNA_CHR$MP24/RNA_CHR$Length*1000
RNA_CHR$RPK_MP48<-RNA_CHR$MP48/RNA_CHR$Length*1000
RNA_CHR$RPK_Y24<-RNA_CHR$Y24/RNA_CHR$Length*1000
RNA_CHR$RPK_Ycongo24<-RNA_CHR$Ycongo24/RNA_CHR$Length*1000
RNA_CHR$RPK_Y48<-RNA_CHR$Y48/RNA_CHR$Length*1000
# Compile these five new columns into a matrix
mydata <- RNA_CHR[, 44:48]
mymatrix <- as.matrix(mydata)
# Scale the data by row (genes) for better visualization
scaled_data <- scale(mydata, center = TRUE, scale = TRUE)
# Create the heatmap with clustering only for columns (conditions)
heatmap(scaled_data,
Rowv = NA, # No clustering for rows (genes)
Colv = TRUE, # Cluster columns (conditions)
col = colorRampPalette(c("green", "white", "red"))(100), # Use diverging colors
margins = c(10, 10),
main = "RNA-seq Heatmap (Column Clustering, Row Scaled)",
cexCol = 1.2, # Adjust column label size
labRow = NA) # Remove row labels
# Calculate the min, mean, and max of the scaled data
scaled_min <- round(min(scaled_data, na.rm = TRUE), digits = 2)
scaled_mean <- round(mean(scaled_data, na.rm = TRUE), digits = 2)
scaled_max <- round(max(scaled_data, na.rm = TRUE), digits = 2)
# Add the legend with the scaled data range
legend("topright",
legend = c(scaled_min, scaled_mean, scaled_max), # Use scaled data range
fill = colorRampPalette(c("green", "white", "red"))(100)[c(1, 50, 100)], # Match heatmap colors
title = "Expression (normalized reads per kb)",
bty = "n")
In this heatmap we can observe that there is a marked global expression
shift as Streptomyces transition from 24h (exponential phase) to 48h
(stationary phase), turning many genes from green (low) to red (high).
And adding the congo Red at 24h induces a transcriptomic change that
resembles late-stage (48h) stationary growth rather than normal 24h
exponential growth.