Lecture 5
Set up
library(ggplot2)
library(dplyr)
library(tximport)
library(DESeq2)
library(vsn)
library(gplots)
library(hexbin)
library(eulerr)
library(GenomicFeatures)
library(rtracklayer)
library(patchwork)
library(tibble)Import data
Annotation
Load
if (!file.exists('gencode.v36.annotation.gff3.gz')) {
download.file('ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_36/gencode.v36.annotation.gff3.gz', 'gencode.v36.annotation.gff3.gz')
}
g <- import.gff3('gencode.v36.annotation.gff3.gz')
tx2gene <- g[g$type == 'transcript'] %>%
{tibble(txid = .$transcript_id, gene = .$gene_id)} %>%
mutate_all(function(x) sub('\\..*', '', x))GFF3 -> GRanges
head(g)## GRanges object with 6 ranges and 23 metadata columns:
## seqnames ranges strand | source type score phase
## <Rle> <IRanges> <Rle> | <factor> <factor> <numeric> <integer>
## [1] chr1 11869-14409 + | HAVANA gene NA <NA>
## [2] chr1 11869-14409 + | HAVANA transcript NA <NA>
## [3] chr1 11869-12227 + | HAVANA exon NA <NA>
## [4] chr1 12613-12721 + | HAVANA exon NA <NA>
## [5] chr1 13221-14409 + | HAVANA exon NA <NA>
## [6] chr1 12010-13670 + | HAVANA transcript NA <NA>
## ID gene_id gene_type
## <character> <character> <character>
## [1] ENSG00000223972.5 ENSG00000223972.5 transcribed_unproces..
## [2] ENST00000456328.2 ENSG00000223972.5 transcribed_unproces..
## [3] exon:ENST00000456328.. ENSG00000223972.5 transcribed_unproces..
## [4] exon:ENST00000456328.. ENSG00000223972.5 transcribed_unproces..
## [5] exon:ENST00000456328.. ENSG00000223972.5 transcribed_unproces..
## [6] ENST00000450305.2 ENSG00000223972.5 transcribed_unproces..
## gene_name level hgnc_id havana_gene
## <character> <character> <character> <character>
## [1] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## [2] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## [3] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## [4] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## [5] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## [6] DDX11L1 2 HGNC:37102 OTTHUMG00000000961.2
## Parent transcript_id transcript_type
## <CharacterList> <character> <character>
## [1] <NA> <NA>
## [2] ENSG00000223972.5 ENST00000456328.2 processed_transcript
## [3] ENST00000456328.2 ENST00000456328.2 processed_transcript
## [4] ENST00000456328.2 ENST00000456328.2 processed_transcript
## [5] ENST00000456328.2 ENST00000456328.2 processed_transcript
## [6] ENSG00000223972.5 ENST00000450305.2 transcribed_unproces..
## transcript_name transcript_support_level tag
## <character> <character> <CharacterList>
## [1] <NA> <NA>
## [2] DDX11L1-202 1 basic
## [3] DDX11L1-202 1 basic
## [4] DDX11L1-202 1 basic
## [5] DDX11L1-202 1 basic
## [6] DDX11L1-201 NA basic
## havana_transcript exon_number exon_id
## <character> <character> <character>
## [1] <NA> <NA> <NA>
## [2] OTTHUMT00000362751.1 <NA> <NA>
## [3] OTTHUMT00000362751.1 1 ENSE00002234944.1
## [4] OTTHUMT00000362751.1 2 ENSE00003582793.1
## [5] OTTHUMT00000362751.1 3 ENSE00002312635.1
## [6] OTTHUMT00000002844.2 <NA> <NA>
## ont protein_id ccdsid
## <CharacterList> <character> <character>
## [1] <NA> <NA>
## [2] <NA> <NA>
## [3] <NA> <NA>
## [4] <NA> <NA>
## [5] <NA> <NA>
## [6] PGO:0000005,PGO:0000019 <NA> <NA>
## -------
## seqinfo: 25 sequences from an unspecified genome; no seqlengthsTranscript to gene mapping
head(tx2gene)## # A tibble: 6 × 2
## txid gene
## <chr> <chr>
## 1 ENST00000456328 ENSG00000223972
## 2 ENST00000450305 ENSG00000223972
## 3 ENST00000488147 ENSG00000227232
## 4 ENST00000619216 ENSG00000278267
## 5 ENST00000473358 ENSG00000243485
## 6 ENST00000469289 ENSG00000243485featureCounts
Load
# read in featureCounts result
fc <- read.table("ex_featurecounts.txt", header = 1)
# rename counts column for simplicity
names(fc)[7] <- 'count'
# take a look at the table
head(fc)## Geneid Chr
## 1 ENSG00000223972.5 chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1
## 2 ENSG00000227232.5 chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1
## 3 ENSG00000278267.1 chr1
## 4 ENSG00000243485.5 chr1;chr1;chr1;chr1;chr1
## 5 ENSG00000284332.1 chr1
## 6 ENSG00000237613.2 chr1;chr1;chr1;chr1;chr1
## Start
## 1 11869;12010;12179;12613;12613;12975;13221;13221;13453
## 2 14404;15005;15796;16607;16858;17233;17606;17915;18268;24738;29534
## 3 17369
## 4 29554;30267;30564;30976;30976
## 5 30366
## 6 34554;35245;35277;35721;35721
## End
## 1 12227;12057;12227;12721;12697;13052;14409;13374;13670
## 2 14501;15038;15947;16765;17055;17368;17742;18061;18366;24891;29570
## 3 17436
## 4 30039;30667;30667;31109;31097
## 5 30503
## 6 35174;35481;35481;36073;36081
## Strand Length count
## 1 +;+;+;+;+;+;+;+;+ 1735 0
## 2 -;-;-;-;-;-;-;-;-;-;- 1351 1
## 3 - 68 0
## 4 +;+;+;+;+ 1021 0
## 5 + 138 0
## 6 -;-;-;-;- 1219 0Distribution
ggplot(fc, aes(x = count)) +
geom_histogram(bins = 200) +
scale_y_log10()
Salmon
Load
# read in raw data
sf <- read.table('quant.sf', header = 1)
Sys.setenv(PATH = paste("/Users/padilr1/opt/anaconda3/envs/r_env/bin", Sys.getenv("PATH"), sep=":"))
# import salmon results
sm <- tximport("quant.sf",
type = "salmon",
tx2gene = tx2gene,
ignoreTxVersion = T)Per-transcript
head(sf)## Name Length EffectiveLength TPM NumReads
## 1 ENST00000456328.2 1657 1406.00 0.000000 0
## 2 ENST00000450305.2 632 381.00 0.000000 0
## 3 ENST00000488147.1 1351 1678.64 4.302518 1
## 4 ENST00000619216.1 68 2.00 0.000000 0
## 5 ENST00000473358.1 712 461.00 0.000000 0
## 6 ENST00000469289.1 535 284.00 0.000000 0Per-gene
str(sm)## List of 4
## $ abundance : num [1:5437, 1] 1.4 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5437] "ENSG00000000457" "ENSG00000000460" "ENSG00000000938" "ENSG00000000971" ...
## .. ..$ : NULL
## $ counts : num [1:5437, 1] 1 0 0 0 0 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5437] "ENSG00000000457" "ENSG00000000460" "ENSG00000000938" "ENSG00000000971" ...
## .. ..$ : NULL
## $ length : num [1:5437, 1] 5141 2180 1471 2309 1434 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:5437] "ENSG00000000457" "ENSG00000000460" "ENSG00000000938" "ENSG00000000971" ...
## .. ..$ : NULL
## $ countsFromAbundance: chr "no"Comparison
Merge
# merge by ensembl gene id
both <- fc %>% mutate(Geneid = sub('\\..*', '', Geneid)) %>%
merge(sm$counts, by.x = 'Geneid', by.y = 'row.names') %>%
dplyr::rename(sm = V1, fc = count)Unexpressed
fit <- lapply(both[c('sm', 'fc')], `==`, 0) %>%
bind_cols() %>%
euler()
plot(fit, quantities = T)
Expressed
both %>%
dplyr::filter(fc != 0 & sm != 0) %>%
ggplot(aes(x = sm, y = fc)) +
geom_point() +
scale_x_log10() +
scale_y_log10() +
labs(x = 'Salmon', y = 'featureCounts')
Create DESeq object
mat <- as.matrix(both[,c('fc','sm')])
cdat <- data.frame(kind = colnames(mat)) %>% `rownames<-`(colnames(mat))
dds <- DESeqDataSetFromMatrix(countData = round(mat),
colData = cdat,
design = ~kind)
dds## class: DESeqDataSet
## dim: 5437 2
## metadata(1): version
## assays(1): counts
## rownames: NULL
## rowData names(0):
## colnames(2): fc sm
## colData names(1): kindFull dataset
Load
# load full dataset
load('lec5.rda')
# take a look at "data"
head(data)## Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3
## ENSG00000000003 4376 3646 2968 6177 3688
## ENSG00000000005 27 29 23 37 13
## ENSG00000000419 1478 1288 884 2368 1340
## ENSG00000000457 508 405 357 934 571
## ENSG00000000460 1394 1165 850 2172 1217
## ENSG00000000938 1 1 0 2 0
## Mov10_oe_1 Mov10_oe_2 Mov10_oe_3
## ENSG00000000003 3343 3114 2079
## ENSG00000000005 25 38 22
## ENSG00000000419 1890 1766 1271
## ENSG00000000457 647 591 354
## ENSG00000000460 1184 1139 673
## ENSG00000000938 0 0 0Variance stabilization
Although differential expression analysis requires raw counts, transformed data generally works better for visualization, clustering, and other uses. Compare a simple log2 transformation versus the regularized log and variance stabilizing transformations.
keep <- rowSums(counts(dds)) >= 10
dds.filt <- dds[keep,]
# variance stabilizing, regularized log, and shift log transformations
lapply(c('vst', 'rlog', 'normTransform'), function(x) {
match.fun(x)(dds.filt) %>%
assay() %>%
meanSdPlot(plot = F) %>%
.$gg +
ggtitle(x)
}) %>% wrap_plots(nrow = 1)
Dimension reduction
vsd <- vst(dds.filt)
plotPCA(vsd)
Hierarchical clustering
There are too many genes to perform hierarchical clustering quickly, so a subset of data will be looked at instead.
# how can you subset the data?
select <- assay(vsd) %>%
apply(1, sd) %>%
order(decreasing = TRUE) %>%
.[1:100]
heatmap(assay(vsd)[select,], margins = c(10, 5))
Spot-check
Pick a gene, does the expression level seem different between conditions?
plotCounts(dds.filt, gene = select[1], intgroup = "condition")
Your turn
Load the featureCount dataset ex2_featureCounts.txt.
Load
fc <- read.table("~/Documents/HGEN_663/extra/lec5/ex2_featureCounts.txt",header=1)Create DESeq object
Remove unnecessary columns. An example of how to generate a metadata
object is indicated below. Within the metadata dataframe, create another
column called condition, grouping together the parental
(PA), NSD1 knockout (NSD1KO) and NSD1/2 double knockouts (NSD12DKO)
samples. Use ~condition in the design parameter when
creating your DESeq object.
df <- fc %>%
dplyr::select(-c(2:6)) %>%
tibble::column_to_rownames("Geneid")
mat <- df %>% as.matrix()
metadata <- data.frame(kind=colnames(mat))
metadata$condition <- "PA"
metadata[3:4,"condition"] <- "NSD1KO"
metadata[5:6,"condition"] <- "NSD12DKO"
metadata <- tibble::column_to_rownames(metadata,"kind")
dds <- DESeqDataSetFromMatrix(countData = mat,
colData = metadata,
design = ~condition)Apply vst transformation
vst <- varianceStabilizingTransformation(dds)PCA
Plot only the top 500 most variable genes.
# subset data
select <- assay(vst) %>%
apply(1, sd) %>%
order(decreasing = TRUE) %>%
.[1:500]
vst_filt <- vst[select,]
#
plotPCA(vst_filt)
Hierarchical clustering
Do the same as above for hierarchical clustering.
heatmap(assay(vst_filt), margins = c(10, 5))
Spotcheck
Google the ensembl id for NSD1 (in humans). Plot the expression level. Does the expression level seem different between the knockout and parental samples?
# NSD1 = ENSG00000165671
# NSD2 = ENSG00000109685
plotCounts(dds,gene = "ENSG00000165671",intgroup = "condition")