RNA_SEQ_TranscriptionFactor_seedling

#library calling
library("dplyr")
library("DESeq2")
library("ggplot2")

#Importing the count data
setwd("C:/Users/user/Desktop/MY Thesis journey(PBT)/RNA seq")
count_data<-read.csv('count_matrix_all_TF.csv',row.names = 1,header=TRUE)
count_data<-count_data[which(rowSums(count_data)>50),]
head(count_data)

#replace with non_zero value for creating deseq object which doesn't except non-Zero value
count_data<-count_data+1
head(count_data)

#importing metadata and filtering for leaf tissue
metadata<-read.csv('metadata_processed.csv',row.names =1,header=TRUE)
#%>%filter(tissue=="leaf")
head(metadata)

#Make DESeq2 object from counts and metadata
ddsMat <- DESeqDataSetFromMatrix(countData = count_data,
                                 colData = metadata,
                                 design = ~treatment)

converting counts to integer mode

ddsMat

class: DESeqDataSet 
dim: 1582 478 
metadata(1): version
assays(1): counts
rownames(1582): Os01g0102400 Os01g0108400 ... Os12g0636200.2 Os12g0636200.3
rowData names(0):
colnames(478): SRR7723492 SRR7723493 ... SRR7723969 SRR7723970
colData names(6): Run Isolation_source ... RNAseqTime_Point PHENOTYPE

# Convert all samples to vsd(data normalization)
vsd <- varianceStabilizingTransformation(ddsMat, blind = FALSE)

-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.

#PCA PLOT
# Plot PCA by column variable
plotPCA(vsd, intgroup = "treatment", ntop = 5) +
  theme(aspect.ratio=1) 

  # remove default ggplot2 theme
  geom_point(size = 1) + 
  # Increase point size
  scale_y_continuous(limits = c(-15, 100)) + 
  # change limits to fix figure dimensions
  ggtitle(label = "Principal Component Analysis (PCA)", 
          subtitle = "Top 500 most variable genes")

Error in `+.gg`:
! Cannot add ggproto objects together. Did you forget to add this object to a ggplot object?
Backtrace:
 1. ggplot2:::`+.gg`(...)

# Find differential expressed genes

ddsMat <- DESeq(ddsMat)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 7 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

ddsMat

class: DESeqDataSet 
dim: 1582 478 
metadata(1): version
assays(6): counts mu ... replaceCounts replaceCooks
rownames(1582): Os01g0102400 Os01g0108400 ... Os12g0636200.2 Os12g0636200.3
rowData names(23): baseMean baseVar ... maxCooks replace
colnames(478): SRR7723492 SRR7723493 ... SRR7723969 SRR7723970
colData names(8): Run Isolation_source ... sizeFactor replaceable

# Get results from testing with FDR adjust pvalues 
results <- results(ddsMat, pAdjustMethod = "fdr", alpha = 0.05)
results

log2 fold change (MLE): treatment stress vs control 
Wald test p-value: treatment stress vs control 
DataFrame with 1582 rows and 6 columns
                baseMean log2FoldChange     lfcSE      stat     pvalue      padj
               <numeric>      <numeric> <numeric> <numeric>  <numeric> <numeric>
Os01g0102400     1.43934      0.1643464  0.160625   1.02317 0.30622921        NA
Os01g0108400     1.49675      0.0782894  0.159634   0.49043 0.62382956        NA
Os01g0108600     3.83023     -0.3451047  0.112826  -3.05872 0.00222282 0.0105697
Os01g0111500     1.61228      0.2442828  0.151661   1.61072 0.10724060 0.2285402
Os01g0128000     1.53146      0.3674500  0.148828   2.46895 0.01355099        NA
...                  ...            ...       ...       ...        ...       ...
Os12g0632600     1.49242     -0.0178436  0.154581 -0.115432   0.908103        NA
Os12g0636200    18.94205      0.0220489  0.129714  0.169980   0.865026  0.930429
Os12g0636200.1  18.94205      0.0220489  0.129714  0.169980   0.865026  0.930429
Os12g0636200.2  18.94205      0.0220489  0.129714  0.169980   0.865026  0.930429
Os12g0636200.3  18.94205      0.0220489  0.129714  0.169980   0.865026  0.930429

# Generate summary of testing.  
summary(results)

out of 1582 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 144, 9.1%
LFC < 0 (down)     : 275, 17%
outliers [1]       : 0, 0%
low counts [2]     : 184, 12%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

resFilt <- results[which(results$padj < 0.05 & abs(results$log2FoldChange) < 1), ]
summary(resFilt)

out of 392 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 130, 33%
LFC < 0 (down)     : 262, 67%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

plotMA(results)

#Volcano plot
res1 = as.data.frame(results)
res1

# add a new column using the mutate function in dplyr
res1 = mutate(res1, sig=ifelse(res1$padj<0.05, "FDR<0.05", "Not Sig"))
res1[which(abs(res1$log2FoldChange)<0.5),'sig'] <- "Not Sig"
res1

p = ggplot(res1, aes(log2FoldChange, -log10(padj))) +
  geom_point(aes(col=sig)) +
  scale_color_manual(values=c("red", "black"))
p

Warning: Removed 184 rows containing missing values (geom_point).

#Gene Annotations
#Plot counts of top gene

topGene <- rownames(res1)[1]
plotCounts(ddsMat, gene = topGene, intgroup= "treatment")