ptc-manuscript-figures/fig_D.Rmd at main · uab-cgds-worthey/ptc-manuscript-figures · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
----
title: "deseq2_subtypes_after-swap.Rmd"
author: "Gurpreet Kaur"
output: html_document
---

```{r}
# a) PCA among 6 subtypes

## DESeq2
dds_sub_after_swap = DESeqDataSetFromMatrix(countData = round(counts), colData = samples_all4deseq_after_swap, design = ~ Gender + Batch +  Phenotype_Subtypes)

keep = rowSums(counts(dds_sub_after_swap)) >= 10
dds_sub_after_swap = dds_sub_after_swap[keep,]

dds_sub_after_swap =  DESeq(dds_sub_after_swap)
resultsNames(dds_sub_after_swap)

## Transformation
vsd_sub_after_swap = vst(dds_sub_after_swap, blind=FALSE)
counts_vst_sub_after_swap = assay(vsd_sub_after_swap)

## Batch correction
library(limma)
counts_vst_sub_after_swap = assay(vsd_sub_after_swap)
mm_sub = model.matrix(~ Gender + Phenotype_Subtypes, colData(vsd_sub_after_swap))
counts_vst_sub_limma_after_swap = limma::removeBatchEffect(counts_vst_sub_after_swap, batch=vsd_sub_after_swap$Batch, design=mm_sub)
assay(vsd_sub_limma_after_swap) = counts_vst_sub_limma_after_swap

## PCA plot
library(DESeq2)
library(dplyr)
pcaData_sub_limma = plotPCA(vsd_sub_limma_after_swap, intgroup=c("Gender", "Batch", "Phenotype_Subtypes"), returnData=TRUE)
percentVar_sub_limma = round(100 * attr(pcaData_sub_limma, "percentVar"))
pcaData_sub_limma = pcaData_sub_limma %>% rename(Sex = Gender)

png("./deseq2_after-swap/results/batch-corr_limma/subtypes/pca_label_batch-corr_subtypes_pid_font.png", width = 4600, height = 3000, res=300)
library(ggplot2)
library(ggrepel)

### Add participant_id
pid = samples_all4deseq_after_swap_pid$Participant_ID[match(samples_all4deseq_after_swap$samples_id, samples_all4deseq_after_swap_pid$samples_id)]

pcaData_sub_limma$Phenotype_Subtypes = factor(pcaData_sub_limma$Phenotype_Subtypes)
pcaData_sub_limma$Batch = factor(pcaData_sub_limma$Batch)
pcaData_sub_limma$Sex = factor(pcaData_sub_limma$Sex)

ggplot(pcaData_sub_limma, aes(PC1, PC2, fill=Phenotype_Subtypes, shape=Batch, color=Sex)) +
  scale_color_manual(values = c("indianred","steelblue1") ) +
  scale_shape_manual(values = c(24, 21, 23)) +
  guides(fill=guide_legend(override.aes=list(shape=21))) +
  guides(color=guide_legend(override.aes=list(shape=21, color=c("indianred","steelblue1")))) +
  scale_fill_brewer(palette="Paired")  +
  guides(
    fill = guide_legend(order = 1),
    shape = guide_legend(order = 2),
    color = guide_legend(order = 3)
  ) +
  geom_label_repel(aes(label = pid), size=5, vjust = "inward", nudge_y=1, max.overlaps = Inf) +
  geom_point(size=2) + xlab(paste0("PC1: ",percentVar_sub_limma[1],"% variance")) + ylab(paste0("PC2: ",percentVar_sub_limma[2],"% variance")) + coord_fixed() + theme_bw(base_size = 20)
dev.off()
```

```{r}
# b) Upset plot for five contrasts
upsetList_5comp_sig = list(
                          res_tum_vs_norm_05 = res_tum_vs_norm_05$gene_name,
                          res_tum_ptc_vs_norm_ptc_05 = res_tum_ptc_vs_norm_ptc_05$gene_name,
                          res_tum_ftc_vs_norm_ftc_05 = res_tum_ftc_vs_norm_ftc_05$gene_name,
                          res_tum_ptc_vs_tum_ftc_05  = res_tum_ptc_vs_tum_ftc_05$gene_name,
                          res_tum_ptc_and_ptcplusthy_vs_tum_ftc_05 = res_tum_ptc_and_ptcplusthy_vs_tum_ftc_05$gene_name)

library(genekitr)
png("./deseq2_after-swap/results/batch-corr_limma/subtypes/tum_ptc_and_ptcplusthy_vs_tum_ftc/upsetgkplot_5comp_sig_color.png", width = 4600, height = 3000, res=300)
plotVenn(upsetList_5comp_sig,
         use_venn = FALSE,
         main_text_size = 15,
         legend_text_size = 8,
         legend_position = 'left',
         color = c("grey70", "#CAB2D6","#FDBF6F", "#B15928", "#6A3D9A")
)
dev.off()
```