BulkRNAseqAnalysis_from_countMatrix/functions.R at main · moritzmadern/BulkRNAseqAnalysis_from_countMatrix · 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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
## functions required for "bulkRNAseq_Analysis.Rmd"


PCA_plot <- function(m, groups, batch, legend_colors, plot_path=NULL, comp1=1, comp2=2){

  # replace NAs with 0
  m[is.na(m)] <- 0

  # berechne PCA
  pca_res <- prcomp(t(m))
  rot_mat <- pca_res$rotation
  res_final <- as.matrix(scale(t(m), center=TRUE, scale=FALSE)) %*% rot_mat
  eigenv <- pca_res$sdev^2
  fraction_var_pca1 <- round(eigenv[comp1]/sum(eigenv),digits=3)
  fraction_var_pca2 <- round(eigenv[comp2]/sum(eigenv),digits=3)

  ## create groups
  groups <- as.factor(groups)
  colors <- legend_colors
  names(colors) <- levels(groups)

  ## ggplot PCR
  if(is.null(batch)){
    df_gg <- as.data.frame(res_final)
    df_gg$samplenames <- samplenames
    df_gg$groups <- groups
    gg <- ggplot(df_gg) +
      geom_point(aes(x=res_final[,comp1], y=res_final[,comp2], col=groups, text=samplenames),size=5) +
      scale_color_manual(values=legend_colors)+
      xlab(paste0("PC",comp1, " (",fraction_var_pca1*100,"%",")")) +
      ylab(paste0("PC", comp2,   " (",fraction_var_pca2*100,"%",")")) +
      theme_bw()
  } else {
    df_gg <- as.data.frame(res_final)
    df_gg$samplenames <- samplenames
    df_gg$groups <- groups
    df_gg$batch <- as.factor(batch)
    gg <- ggplot(df_gg) +
      geom_point(aes(x=res_final[,comp1], y=res_final[,comp2], col=groups, text=samplenames, shape=batch),size=5) +
      scale_color_manual(values=legend_colors)+
      xlab(paste0("PC",comp1, " (",fraction_var_pca1*100,"%",")")) +
      ylab(paste0("PC", comp2,   " (",fraction_var_pca2*100,"%",")")) +
      theme_bw()
  }

  # save plot
  ggsave(plot=gg, filename=plot_path, width = 6, height = 4)

  # print plot
  ggplotly(gg)

}


### heatmap plot ###
heatmap_plot <- function(m, groups, legend_colors, sample_names, type="normal", dendrogram="column", labrow="", bool_rowv=TRUE, bool_colv = TRUE){

  # create groups
  names(colors) <- levels(groups)

  # replaces NAs with 0
  m[is.na(m)] <- 0
  colnames(m) <- sample_names

  # should rows be reordered
  if (bool_rowv){
    rowv <- as.dendrogram(hclust(dist(m)))
  } else {
    rowv <- FALSE
  }

  # should columns be reordered
  if (is.logical(bool_colv)){
    if (bool_colv){
      colv <- as.dendrogram(hclust(dist(t(m))))
    } else {
      colv <- FALSE
    }
  } else {
    colv <- bool_colv
  }

  # specify colors
  if(is.null(legend_colors)){
    sidecolors <- rep("white", times=ncol(m))
  } else{
    sidecolors <- legend_colors[groups]
  }

  # create color palette
  colors_heatmap <- rev(brewer.pal(11, "RdBu"))
  colors_heatmap[6] <- "#fffec8"
  heatmap_pal <- colorRampPalette(colors_heatmap)
  rdbu_colors = heatmap_pal(20)[2:19]
  heatmap_pal <- colorRampPalette(rdbu_colors)

  # plot heatmap
  colors_pal <- colorRampPalette(pals::parula(40))
  par(mfrow=c(1,1))
  par(xpd=TRUE)
  if (type == "normal"){
    heatmap.2(m,
              Rowv = rowv,
              Colv=colv,
              margins=c(8,8), cexCol = 1,labRow=labrow,col=heatmap_pal(50), ColSideColors = sidecolors, symkey = F,
              cex.lab=1.5, scale="none", trace="none", dendrogram=dendrogram)
  }
  if (type == "centered"){
    min_m <- min(m, na.rm=TRUE)
    max_m <- max(m, na.rm=TRUE)
    heatmap.2(m,
              Rowv = rowv,
              Colv=colv,
              labRow=labrow, margins=c(8,8), ColSideColors = sidecolors, trace="none",col=heatmap_pal(50),
              breaks = seq(from=-2,to=2, length.out=51),
              symkey = F,
              dendrogram=dendrogram)
  }
  if (type == "standardized"){
    min_m <- min(m, na.rm=TRUE)
    max_m <- max(m, na.rm=TRUE)
    heatmap.2(m,
              Rowv = rowv,
              Colv= colv,
              labRow=labrow, margins=c(8,8), ColSideColors = sidecolors, trace="none",col=heatmap_pal(50), symkey = F,
              breaks = seq(from=-2,to=2, length.out=51),
              dendrogram=dendrogram)
  }

  # add legend
  par(xpd=TRUE)
}