DESeq2 Analysis - Dataset 1

#calling the libraries

library(DESeq2)
library(ggplot2)

#calling in the count matrix
counts_D1 <- read.delim("counts_D1", header = TRUE, row.names = 1)

#understanding the distribution of counts in the data set

total_counts_per_gene_D1 <- rowSums(counts_D1)

hist(total_counts_per_gene_D1, 
     breaks=100, 
     main="Distribution of Total Counts per Gene", 
     xlab="Total counts across all samples",
     ylab="Number of genes")

summary(total_counts_per_gene_D1)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0       0       6    3657     191 3392661 

#ignoring low expressed gene to remove noise & get better results

counts_D1 <- counts_D1[which(rowSums(counts_D1) > 50), ]

#calling in the metadata
metadata_D1 <- read.csv("metadata_D1.csv", header = TRUE)

#rounding of the decimals, DeSeq2 cannot handle them -> creating the DeSeq2 object

counts_D1 <- round(counts_D1)
dds_D1 <- DESeqDataSetFromMatrix(countData = counts_D1,
                             colData = metadata_D1,        
                             design = ~Treatment)

dds_D1

class: DESeqDataSet 
dim: 25946 115 
metadata(1): version
assays(1): counts
rownames(25946): ENSG00000142611 ENSG00000157911 ... ENSG00000307722
  ENSG00000302039
rowData names(0):
colnames(115): SRR26115245 SRR26115246 ... SRR26115358 SRR26115359
colData names(2): Run Treatment

colData(dds_D1)

DataFrame with 115 rows and 2 columns
                    Run      Treatment
            <character>       <factor>
SRR26115245 SRR26115245 Osimertinib 6w
SRR26115246 SRR26115246 Osimertinib 6w
SRR26115247 SRR26115247 Osimertinib 6w
SRR26115248 SRR26115248 Osimertinib 6w
SRR26115249 SRR26115249 Osimertinib 6w
...                 ...            ...
SRR26115355 SRR26115355      Untreated
SRR26115356 SRR26115356      Untreated
SRR26115357 SRR26115357      Untreated
SRR26115358 SRR26115358      Untreated
SRR26115359 SRR26115359      Untreated

#THIS IS THE CORE ANALYSIS
dds_D1 <- DESeq(dds_D1) 

#Variance stabilizing transformation for visualization
vsdata_D1 <- vst(dds_D1, blind = TRUE) 

#plotting the PCA, got good results
plotPCA(vsdata_D1, intgroup = "Treatment") 

#results name
resultsNames(dds_D1) 

[1] "Intercept"                             "Treatment_Untreated_vs_Osimertinib.6w"

#Extracts the final DE genes
res_D1 <- results(dds_D1) 

#Summarizes results
summary(res_D1) 

out of 25946 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 6600, 25%
LFC < 0 (down)     : 8491, 33%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

#plotting the dispersion plot
plotDispEsts(dds_D1)

What the dispersion plot means for your analysis:

• DESeq2’s statistical model fits your data well
• Your p-values and fold changes are reliable
• No major technical issues or batch effects
• Good quality RNA-seq data

#getting the results name & making contrast
resultsNames(dds_D1)

[1] "Intercept"                             "Treatment_Untreated_vs_Osimertinib.6w"

res_D1 <- results(dds_D1, contrast = c("Treatment", "Osimertinib 6w", "Untreated"))
res_D1

log2 fold change (MLE): Treatment Osimertinib 6w vs Untreated 
Wald test p-value: Treatment Osimertinib 6w vs Untreated 
DataFrame with 25946 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat      pvalue        padj
                <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
ENSG00000142611  13.36936     -3.0607931 0.3574999 -8.561662 1.11252e-17 1.41082e-16
ENSG00000157911  30.48609      0.0808536 0.0623283  1.297222 1.94555e-01 2.78015e-01
ENSG00000142655  34.30945     -0.1170707 0.0508446 -2.302521 2.13058e-02 4.18186e-02
ENSG00000149527   9.25249     -0.3146144 0.3724180 -0.844788 3.98229e-01 4.97757e-01
ENSG00000171621  30.82709     -0.4951862 0.0967244 -5.119556 3.06255e-07 1.42993e-06
...                   ...            ...       ...       ...         ...         ...
ENSG00000298139  3.005739       1.161385  0.192194   6.04277 1.51488e-09 9.30636e-09
ENSG00000298153  0.418607       1.764653  0.399907   4.41266 1.02107e-05 3.85419e-05
ENSG00000298169 14.029668       0.606413  0.125423   4.83495 1.33176e-06 5.67202e-06
ENSG00000307722  0.623014      -1.751171  0.320955  -5.45613 4.86630e-08 2.50468e-07
ENSG00000302039  3.781707       0.249161  0.168233   1.48105 1.38594e-01 2.09750e-01

summary(res_D1)

out of 25946 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 8491, 33%
LFC < 0 (down)     : 6600, 25%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

library(apeglm)

#applying LFC shrinkage
res_shrunk_D1 <- lfcShrink(dds_D1, coef=2, type="apeglm")  # coef=2 is treatment effect

#significant genes from SHRUNKEN results
sig_overall_D1 <- res_shrunk_D1[which(res_shrunk_D1$padj < 0.05 & abs(res_shrunk_D1$log2FoldChange) > 1), ]

#the top & down regulated gene lists
top_up_D1 <- head(sig_overall_D1[sig_overall_D1$log2FoldChange > 1,][order(sig_overall_D1[sig_overall_D1$log2FoldChange > 1,]$padj),], 10)
top_down_D1 <- head(sig_overall_D1[sig_overall_D1$log2FoldChange < -1,][order(sig_overall_D1[sig_overall_D1$log2FoldChange < -1,]$padj),], 10)

#Export the results
write.csv(as.data.frame(res_shrunk_D1), "D1_full_results.csv")
write.csv(sig_overall_D1, "significant_genes_D1.csv")
write.csv(counts(dds_D1, normalized=TRUE), "normalized_counts_D1.csv")
write.csv(assay(vsdata_D1), "vst_data_D1.csv")

Visualization Plots

Pathway Analysis

# Gene Ontology Analysis
BiocManager::install("clusterProfiler")

Update all/some/none? [a/s/n]: 

a

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/cli_3.6.5.zip'
Content type 'application/zip' length 1402276 bytes (1.3 MB)
downloaded 1.3 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/jsonlite_2.0.0.zip'
Content type 'application/zip' length 1109640 bytes (1.1 MB)
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/rlang_1.1.6.zip'
Content type 'application/zip' length 1629108 bytes (1.6 MB)
downloaded 1.6 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/xfun_0.52.zip'
Content type 'application/zip' length 607535 bytes (593 KB)
downloaded 593 KB

package ‘cli’ successfully unpacked and MD5 sums checked
package ‘jsonlite’ successfully unpacked and MD5 sums checked
package ‘rlang’ successfully unpacked and MD5 sums checked
package ‘xfun’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\aasth\AppData\Local\Temp\Rtmpm6j6Xe\downloaded_packages

library(clusterProfiler)
BiocManager::install("org.Hs.eg.db")

Update all/some/none? [a/s/n]: 

a

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/cli_3.6.5.zip'
Content type 'application/zip' length 1402276 bytes (1.3 MB)
downloaded 1.3 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/jsonlite_2.0.0.zip'
Content type 'application/zip' length 1109640 bytes (1.1 MB)
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/rlang_1.1.6.zip'
Content type 'application/zip' length 1629108 bytes (1.6 MB)
downloaded 1.6 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/xfun_0.52.zip'
Content type 'application/zip' length 607535 bytes (593 KB)
downloaded 593 KB

package ‘cli’ successfully unpacked and MD5 sums checked
package ‘jsonlite’ successfully unpacked and MD5 sums checked
package ‘rlang’ successfully unpacked and MD5 sums checked
package ‘xfun’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\aasth\AppData\Local\Temp\Rtmpm6j6Xe\downloaded_packages

library(org.Hs.eg.db)

# Get significant upregulated genes (potential resistance mechanisms)
up_genes_D1 <- rownames(res_shrunk_D1[res_shrunk_D1$log2FoldChange > 1 & 
                               res_shrunk_D1$padj < 0.05 & 
                               !is.na(res_shrunk_D1$padj), ])

# Get significant downregulated genes (potential targets)
down_genes_D1 <- rownames(res_shrunk_D1[res_shrunk_D1$log2FoldChange < -1 & 
                                 res_shrunk_D1$padj < 0.05 & 
                                 !is.na(res_shrunk_D1$padj), ])

# GO enrichment for upregulated genes
go_up_D1 <- enrichGO(gene = up_genes_D1,
                 OrgDb = org.Hs.eg.db,
                 keyType = "ENSEMBL",
                 ont = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05)

dotplot(go_up_D1, showCategory = 15, title = "Upregulated Pathways")

# GO enrichment for downregulated genes  
go_down_D1 <- enrichGO(gene = down_genes_D1,
                   OrgDb = org.Hs.eg.db,
                   keyType = "ENSEMBL", 
                   ont = "BP",
                   pAdjustMethod = "BH",
                   pvalueCutoff = 0.05)

dotplot(go_down_D1, showCategory = 15, title = "Downregulated Pathways")

GSEA Analysis (Gene Set Enrichment Analysis)

# Prepare ranked gene list for GSEA
BiocManager::install("fgsea")

Update all/some/none? [a/s/n]: 

a

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/cli_3.6.5.zip'
Content type 'application/zip' length 1402276 bytes (1.3 MB)
downloaded 1.3 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/jsonlite_2.0.0.zip'
Content type 'application/zip' length 1109640 bytes (1.1 MB)
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/rlang_1.1.6.zip'
Content type 'application/zip' length 1629108 bytes (1.6 MB)
downloaded 1.6 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/xfun_0.52.zip'
Content type 'application/zip' length 607535 bytes (593 KB)
downloaded 593 KB

package ‘cli’ successfully unpacked and MD5 sums checked
package ‘jsonlite’ successfully unpacked and MD5 sums checked
package ‘rlang’ successfully unpacked and MD5 sums checked
package ‘xfun’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\aasth\AppData\Local\Temp\Rtmpm6j6Xe\downloaded_packages

BiocManager::install("msigdbr")

Update all/some/none? [a/s/n]: 

a

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/cli_3.6.5.zip'
Content type 'application/zip' length 1402276 bytes (1.3 MB)
downloaded 1.3 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/jsonlite_2.0.0.zip'
Content type 'application/zip' length 1109640 bytes (1.1 MB)
downloaded 1.1 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/rlang_1.1.6.zip'
Content type 'application/zip' length 1629108 bytes (1.6 MB)
downloaded 1.6 MB

trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.4/xfun_0.52.zip'
Content type 'application/zip' length 607535 bytes (593 KB)
downloaded 593 KB

package ‘cli’ successfully unpacked and MD5 sums checked
package ‘jsonlite’ successfully unpacked and MD5 sums checked
package ‘rlang’ successfully unpacked and MD5 sums checked
package ‘xfun’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\aasth\AppData\Local\Temp\Rtmpm6j6Xe\downloaded_packages

library(fgsea)
library(msigdbr)

# Create ranked gene list (all genes, ranked by stat)
gene_ranks_D1 <- res_shrunk_D1$log2FoldChange
names(gene_ranks_D1) <- rownames(res_shrunk_D1)
gene_ranks_D1 <- gene_ranks_D1[!is.na(gene_ranks_D1)]
gene_ranks_D1 <- sort(gene_ranks_D1, decreasing = TRUE)

# Get Hallmark pathways
hallmark_sets_D1 <- msigdbr(species = "Homo sapiens", category = "H")
hallmark_list_D1 <- split(hallmark_sets_D1$ensembl_gene, hallmark_sets_D1$gs_name)

# Run GSEA
fgsea_results_D1 <- fgsea(pathways = hallmark_list_D1,
                      stats = gene_ranks_D1,
                      minSize = 15,
                      maxSize = 500,
                      nperm = 1000)

# View results
head(fgsea_results_D1[order(pval), ], 15)

# Plot top pathways
plotEnrichment(hallmark_list_D1[["HALLMARK_APOPTOSIS"]], gene_ranks_D1) + 
  labs(title = "Apoptosis Pathway")