PCA_script
Francisco Carranza
2024-03-19
Franco / DEseq2 PCA
Overview
This a walk-through of the PCA plot script used for Met Cancer project in the Hertel and Kaiser lab. There is a built in PCA command in the DEseq2 R package. For a Deseq2 tutorial/walkthrough see https://altsplicer.github.io/DEseq2_Script/DESEQ2_met.html. You can also see this walk-through via this link. # PCA plots are used as QC check in your RNA-seq analysis. To explore the similarity of our samples and identify outliers. # Highly recommended to do prior to your differential expression analysis. # For a more indepth view of QC checks such as PCA plots and unsupervised clustering see the following link.
Install and Load DEseq2 and ggplot2 in R studio
# installion of DESEQ2 via Biocmanager and ggplot2
#BiocManager::install("DESeq2")
#install.packages("ggplot2")
# Call required packages
library(DESeq2)## Warning: package 'DESeq2' was built under R version 4.3.3## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following object is masked from 'package:grDevices':
##
## windows## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedianslibrary(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3Locate your reads count file
The input file will be a the output text file from the STAR aligner, geneCounts or a feature counts text file.
# location of read counts text file
# infile = "./example.txt"
infile = "./com_0v120v720counts.txt"Setup the experimental design
This will set up the experimental design that will be used to run the DEseq2 analysis. In this case we have 3 reps per condition and 8 conditions. 2 cells lines under 4 media conditions.
# design for 3 reps, 6 conditions
design = c( "3", "3", "3", "3", "3", "3")
# design for 2 reps 4 conditions
# design = c( "2", "2", "2", "2")
# call the design and transform them into a integer character.
reps1 = as.integer(design[1])
reps2 = as.integer(design[2])
reps3 = as.integer(design[3])
reps4 = as.integer(design[4])
reps5 = as.integer(design[5])
reps6 = as.integer(design[6])
# Set up the conditions based on the experimental setup.
# cond1 = rep("cond1", reps1)
# cond2 = rep("cond2", reps2)
cond1 = rep("MB_000", reps1)
cond2 = rep("R8_000", reps2)
cond3 = rep("MB_120", reps3)
cond4 = rep("R8_120", reps4)
cond5 = rep("MB_720", reps5)
cond6 = rep("R8_720", reps6)Read the data from the standard input.
counts = read.table(infile, header=TRUE, sep="\t", row.names=1 )
# Optional, head the counts to view the formatting and adjust accordingly
#head(counts, 5)Edit the counts matrix
The count data so the only the read counts and gene name remains.
# Assume the last columns are the count matrix.
#idx = ncol(counts) - (reps1 + reps2)
idx = ncol(counts) - (reps1 + reps2 + reps3 + reps4 + reps5 + reps6)
# Cut out the valid columns.
# counts = counts[-c(1:idx)]
counts = counts[-c(1:idx)]
# Some tools generate the estimated counts as real numbers
# DESeq 2 allows only integers. We need to convert real numbers to rounded integers.
numeric_idx = sapply(counts, mode) == 'numeric'
counts[numeric_idx] = round(counts[numeric_idx], 0)
head(counts, 5)## MB468_000_01 MB468_000_02 MB468_000_03 R8_000_01 R8_000_02 R8_000_03
## C1orf141 0 0 0 0 0 0
## PKP1 48 42 28 79 68 37
## RERE 3487 3719 2234 4710 3530 2740
## CSMD2 0 1 2 0 0 0
## HIVEP3 598 638 350 302 254 193
## MB468_120_01 MB468_120_02 MB468_120_03 R8_120_01 R8_120_02 R8_120_03
## C1orf141 0 0 0 0 0 0
## PKP1 234 170 272 91 175 47
## RERE 2525 1645 2435 3485 5960 2563
## CSMD2 0 0 0 0 0 0
## HIVEP3 447 241 446 314 460 216
## MB468_720_01 MB468_720_02 MB468_720_03 R8_720_01 R8_720_02 R8_720_03
## C1orf141 0 0 0 0 0 0
## PKP1 646 398 562 41 33 25
## RERE 4331 2484 4481 4023 2894 3657
## CSMD2 0 0 0 0 0 0
## HIVEP3 1082 523 1011 381 324 331Build the dataset names and conditions
samples = names(counts)
condition = factor(c(cond1, cond2, cond3, cond4, cond5, cond6))
colData = data.frame(samples=samples, condition=condition)
# You can view the dataset by calling the following
samples## [1] "MB468_000_01" "MB468_000_02" "MB468_000_03" "R8_000_01" "R8_000_02"
## [6] "R8_000_03" "MB468_120_01" "MB468_120_02" "MB468_120_03" "R8_120_01"
## [11] "R8_120_02" "R8_120_03" "MB468_720_01" "MB468_720_02" "MB468_720_03"
## [16] "R8_720_01" "R8_720_02" "R8_720_03"condition## [1] MB_000 MB_000 MB_000 R8_000 R8_000 R8_000 MB_120 MB_120 MB_120 R8_120
## [11] R8_120 R8_120 MB_720 MB_720 MB_720 R8_720 R8_720 R8_720
## Levels: MB_000 MB_120 MB_720 R8_000 R8_120 R8_720colData## samples condition
## 1 MB468_000_01 MB_000
## 2 MB468_000_02 MB_000
## 3 MB468_000_03 MB_000
## 4 R8_000_01 R8_000
## 5 R8_000_02 R8_000
## 6 R8_000_03 R8_000
## 7 MB468_120_01 MB_120
## 8 MB468_120_02 MB_120
## 9 MB468_120_03 MB_120
## 10 R8_120_01 R8_120
## 11 R8_120_02 R8_120
## 12 R8_120_03 R8_120
## 13 MB468_720_01 MB_720
## 14 MB468_720_02 MB_720
## 15 MB468_720_03 MB_720
## 16 R8_720_01 R8_720
## 17 R8_720_02 R8_720
## 18 R8_720_03 R8_720Create DESEq2 dataset.
#dds = DESeqDataSetFromMatrix(countData=counts, colData=colData, design = ~condition)
dds = DESeqDataSetFromMatrix(countData=counts, colData=colData, design = ~condition)## converting counts to integer modeRun deseq2.
dds = DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingTransform counts for PCA
vst will transform your counts for data visualization.
# object refers to a DESeqDataSet or a matrix of counts.
# blind see documentation for an explanation.
#vst(object, blind=TRUE)
vsdata <- vst(dds, blind=TRUE) Plot a simple PCA plot
plotPCA(vsdata, intgroup=c("condition"))## using ntop=500 top features by variance
# Further Customization
Run the previous command but instead of plotting a PCA plot your output will return a data frame and store it in “pcaData”.
# pcaData <- plotPCA(vsdata, intgroup=c("condition", "type"), returnData=TRUE)
pcaData <- plotPCA(vsdata, intgroup=c("condition"), returnData=TRUE)## using ntop=500 top features by variance#retrieve the percent variance for PC1 and PC2 in the PCA plot
percentVar <- round(100 * attr(pcaData, "percentVar"))Plot the PCA plot using ggplot and the “pcaData” dataframe made using the plotPCA command.
# Use ggplot2 options to change the plot colors, font size and many other plot options.
ggplot(pcaData, aes(PC1, PC2, color=condition)) +
geom_point(size=10,alpha=0.5) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
coord_fixed() + scale_color_manual(values = c("#1565C0","#2196F3","#64B5F6","#641E16","#C0392B","#E6B0AA")) +
theme(legend.title=element_blank()) + theme(axis.title=element_text(size=14,face="bold")) + theme(axis.text.x=element_text(size=12))+
theme(axis.text.y=element_text(size=12)) + theme(legend.text=element_text(size=10))
Acknowledgments
- Page header logo is adapted from images: [DESeq2-Michael-Love.png]