deseq_dataset =DESeqDataSetFromMatrix(countData = counts_filtered,colData = design_matrix, # colData means data about the columns of the count matrixdesign =~Genotype # basically: expression in function of the genotype) %>%DESeq()
There are two basic ways we can check for sample homogeneity, meaning how similar the replicates of a sample are to each other: principal component analysis (PCA) or clustering analysis.
To perform either of them, we need to use the normalized data, not the read counts.
Normalize counts
You can either use rlogTransformation or varianceStabilizingTransformation for PCA, as they return an object that includes the metadata. The first takes about 5x longer than the second, but results are similar. Check the help page for details on the computation.
To see the transformed counts, you need to run assay on the result.
rlogTransformation
rlog_transformed =rlogTransformation(deseq_dataset)assay(rlog_transformed) %>% head
DESeq2 offers a function for PCA that takes the result from either of the two functions and returns a ggplot object.
Because of this, we can add things to the plot as if it was a ggplot.
plotPCA(rlog_transformed, intgroup ="Genotype", # color samples byntop =nrow(counts_filtered)) #use how many genes for the computation
PCA
DESeq2 offers a function for PCA that takes the result from either of the two functions and returns a ggplot object.
Because of this, we can add things to the plot as if it was a ggplot.
plotPCA(rlog_transformed, intgroup ="Genotype", # color samples byntop =nrow(counts_filtered)) +#use how many genes for the computationlabs(color ="") +scale_color_manual(values =c("Control"="black","Mlp37347"="red","Mlp124499"="blue") )
Clustering
Hierarchical clustering is a graphical way to show the pairwise similarity between samples.
We just need to edit the rlog_transformed or the VST_data so that the samples are shown in the rows and the genes in the columns.
distance_data = rlog_transformed %>%assay() %>% t %>%dist() # default is euclideancluster_hc =hclust(distance_data)plot(cluster_hc)
Add a heatmap
You can show the dendrogram and the distance matrix in a heatmap
If you want to plot it with ggplot to have more control, you need to do some transformations
# Pivot longer the distance matrix and add the row nameslong_distance <- distance_data %>%as.matrix() %>%as_tibble() %>%# rownames of the distance matrix are in the attr(ibute) Labelsmutate(Rep1 =attr(distance_data, which ="Labels")) %>%pivot_longer(cols =!Rep1,names_to ="Rep2", values_to ="distance")
Improved heatmap + dendrogram: prepare data
If you want to plot it with ggplot to have more control, you need to do some transformations
# Pivot longer the distance matrix and add the row nameslong_distance <- distance_data %>%as.matrix() %>%as_tibble() %>%# rownames of the distance matrix are in the attr(ibute) Labelsmutate(Rep1 =attr(distance_data, which ="Labels")) %>%pivot_longer(cols =!Rep1,names_to ="Rep2", values_to ="distance")# Add a column with the relative distance of the samples(max_distance =max(long_distance$distance))
If you want to plot it with ggplot to have more control, you need to do some transformations
# Pivot longer the distance matrix and add the row nameslong_distance <- distance_data %>%as.matrix() %>%as_tibble() %>%# rownames of the distance matrix are in the attr(ibute) Labelsmutate(Rep1 =attr(distance_data, which ="Labels")) %>%pivot_longer(cols =!Rep1,names_to ="Rep2", values_to ="distance")# Add a column with the relative distance of the samples(max_distance =max(long_distance$distance))
[1] 106.0121
long_distance = long_distance %>%mutate(Rel_distance = (distance / max_distance) *100)# Sort the sample names according to the dendrogramlong_distance$Rep1 =factor(long_distance$Rep1, levels = cluster_hc$labels[cluster_hc$order] )long_distance$Rep2 =factor(long_distance$Rep2, levels = cluster_hc$labels[cluster_hc$order] )
Improved heatmap + dendrogram: plots
Tree should not have the tip labels, cause the labels will be the axis of the heatmap
Normally, we want to plot the deregulated genes, find patterns among them, etc. The first thing we can plot is a volcano plot for each contrast checked.
It is simply the fold change in the x-axis and the -log10p-adj in the y-axis.
Note that there are cases where the p-value, and thus the adjusted p-value, are NA, so we will create a new variable that replaces NA with 1 (those genes are not deregulated)
Now we also want to sort the genes so that the ones with similar deregulation are grouped together. We can use a cluster analysis for that.
But, we need to use a wide table, not a long one, to compute the distance matrix. We will also filter out the genes that are not deregulated, so that the analysis don’t consider non-significant differences.
values_fill = 0: after filtering, there will be genes with fold change levels present for only one contrast because they are not deregulated in the other condition. Since in the wide format, we need values for both contrasts, we replace the empty cell with 0.
Heatmap of deregulated genes - gene clustering
Now we can use wide_results in the dist function?
Yes, after we remove the Gene_name column (dist only accepts numbers)
Whether you plot the genes in the x or in the y axis depends on how readable the plot is in the end (how many genes are deregulated, how many contrasts you have, etc)
results_DESeq2_filtered %>%ggplot(aes(x = orderedGenes, y = shortContrast, fill =ifelse(abs(log2FoldChange) >2& padj <0.01, log2FoldChange, 0) )) +geom_tile()
Heatmap of deregulated genes - plot
Whether you plot the genes in the x or in the y axis depends on how readable the plot is in the end (how many genes are deregulated, how many contrasts you have, etc)
results_DESeq2_filtered %>%ggplot(aes(x = orderedGenes, y = shortContrast, fill =ifelse(abs(log2FoldChange) >2& padj <0.01, log2FoldChange, 0) )) +geom_tile() +scale_fill_gradient2(name ="Fold change",mid ="white", low ="midnightblue", high ="red")
Heatmap of deregulated genes - plot
Whether you plot the genes in the x or in the y axis depends on how readable the plot is in the end (how many genes are deregulated, how many contrasts you have, etc)
results_DESeq2_filtered %>%ggplot(aes(x = orderedGenes, y = shortContrast, fill =ifelse(abs(log2FoldChange) >2& padj <0.01, log2FoldChange, 0) )) +geom_tile() +scale_fill_gradient2(name ="Fold change",mid ="white", low ="midnightblue", high ="red") +theme(axis.title =element_blank(),axis.text.x =element_blank(),legend.position ="bottom",axis.ticks =element_blank())
