mouse_first_look.qmd

---
title: "Ein erster Blick auf die Evo-Devo-Maus-Daten"
format: html
---

Hier verwenden wir Daten aus folgendem Paper:

M. Cardoso-Moreia et al.: *Gene expression across mammalian organ development*. [Nature 571:505 (2019)](https://doi.org/10.1038/s41586-019-1338-5).

Spezifisch arbeiten wir mit der Expressions-Matrix, die aus den Maus-Proben gebildet wurde. Zusammen mit den Rohdaten haben die Autoren diese Matrix als Datei `Mouse.CPM.txt` verfügbar gemacht. Allerdings wurden in dieser Datei die rohen Read-Zahlen bereits in CPM-Werte umgewandelt. Da wir am Anfang beginnen wollen, habe ich diese Normalisierung rückgängig gemacht (siehe [hier](mouse_unnormalize.html) für Details); die Matrix, die ich so erhalten habe, ist hier [hier](Downloads/evodevo_mouse_counts.tsv.gz) verfügbar).

Wir laden diese Matrix

```{r}
suppressPackageStartupMessages( library(tidyverse) )

read_tsv( "Downloads/evodevo_mouse_counts.tsv.gz" ) -> mouse_counts
```

Hier sind die ersten Spalten und Zeilen der Matrix:

```{r}
mouse_counts[ 1:5, 1:5 ]
```

### Probentabelle

Die Tabelle enthält Daten von 316 Proben. Wir verschaffen uns zunächst einen Überblick über diese, indem wir die Spaltennamen interpretieren. Jeder Spaltennamen besteht aus drei Teilen, die durch Unterstriche getrennt sind: dem Gewebe (bei der ersten Spalte: `Brain`), dem Zeitpunkt (`e10.5`) und der Nummer des Replikats (d.h., des Embryos oder Tieres; hier `1`). 

Wir erzeugen eine Tabelle mit den Spaltennamen, die wir als Probennamen ("sample names") verwenden werden. Den Namen der ersten Spalte verwenden wir nicht, da das die Gen-Namen sind.

```{r}
tibble( sample = colnames(mouse_counts)[-1] ) %>% head()
```

Mit dem Tidyverse-Verb `separate` können wir die Inhalte der `sample`-Spalte entlang der Unterstriche aufteilen und in drei neue Spalten schreiben:

```{r}
tibble( sample = colnames(mouse_counts)[-1] ) %>%
separate( sample, c( "tissue", "timepoint", "replicate" ), sep="_", remove=FALSE ) %>%
mutate( timepoint = fct_inorder( timepoint ) ) -> sample_table  

head( sample_table )
```

Hier haben wir `timepoint` zu einem Faktor gemacht, um sicherzustellen, dass die Reihenfolge der Zeitpunkte später beim Plotten nicht verändert wird. Wir haben nun folgende Zeitpunkte:

```{r}
sample_table %>% pull(timepoint) %>% levels()
```

Hier bedeutet "e10.5", dass die Probe von einem Embryo stammt, 10.5 Tage nach der Befruchtung, und "P3", dass sie von einem Tier stammt, 3 Tage nach der Geburt. (Die Funktion `levels` gibt die verschiedenen Werte eine Faktors in der festgelegten Reihenfolge aus.)

Wir haben folgende Gewebearten:

```{r}
sample_table %>% pull( tissue ) %>% unique()
```
(Die Funktion `unique` hier kürzt einen vektor so, dass jeder Wert nur noch einmal vorkommt.)

Der folgende Code zählt, wieviel Proben wir für jede Kombination aus Gewebe und Zeitpunkt haben:

```{r}
sample_table %>%
group_by( tissue, timepoint ) %>%
summarise( n=n() ) %>%
pivot_wider( names_from = "tissue", values_from = "n", values_fill=0 )
```

Erkennen Sie, wie der Code funktioniert? Sehen Sie sich dazu an, wie die Ausgabe aussieht ohne die letzte Zeile. // Für Fortgeschrittene: Was macht wohl das `values_fill`? Was geschieht, wenn man es weglässt, und warum geschieht das?


### Hochexprimierte Gene; Gen-Namen

Nun sollten wir einen ersten Blick auf unsere Count-Matrix werfen:

```{r}
mouse_counts[1:5,1:5]
```

Welche Gene sind wohl in der Leber besonders hoch exprimiert? Wir nehmen uns eine Probe heraus, z.B. die erste Probe mit Leber neugeborener Mäuse, also `Liver_P0_1` und sortieren die Tabelle absteigend ("descending")

```{r}
mouse_counts %>%
select( gene_id, Liver_P0_1 ) %>%
arrange( desc( Liver_P0_1 ) )
```

Mit Google lässt sich schnell ermittlen, dass die Ensembl-ID `ENSMUSG00000029368`, die nun ganz oben steht, für das Gen *Alb* steht, das für Albumin kodiert, das Protein, das im Blutserum mit der höchsten Konzentration vorkommt. Dass die Leberzellen sehr viel mRNA für dieses Gen brauchen, um dieses wichtige Protein in der erfordeelichen großen Menge zu erzeugen, leuchtet ein.

Was sind die anderen Gene? Wir laden die Datei mit der Zuordnung von Ensembl-Gen-IDs zu Gen-Symbolen , die wir bereits das [letzte Mal](alignment.html#gen-ids) erstellt haben. Wir benennen aber die Spalten um, um sie etwas kürzer zu machen:

```{r}
read_tsv( "Downloads/Ensembl_102_GRCm38.p6_Gene_names.tsv" ) %>%
rename( gene_id = "Gene stable ID", gene_name = "Gene name", gene_descr = "Gene description" )-> gene_names

head( gene_names )
```

Wir fügen diese Informationen unserem Ergebnis an, indem wir den Code von eben durch ein `left_join` ergänzen:

```{r}
mouse_counts %>%
select( gene_id, Liver_P0_1 ) %>%
arrange( desc( Liver_P0_1 ) ) %>%
left_join( gene_names ) %>%
head( 10 )
```
Auch die anderen hochexprimierten Gene passen zur Aufgaben der Leber.

Werfen wir zum Vergleich noch einen Blick auf die Gene, die in der Niere am höchsten exprimiert sind. Wir verwenden wieder das erste Replikat von Zeitpunkt P0:

```{r}
mouse_counts %>%
select( gene_id, Kidney_P0_1 ) %>%
arrange( desc( Kidney_P0_1 ) ) %>%
left_join( gene_names ) %>%
head( 10 )
```

Dieses Mal finden wir hauptsächlich basale Transkriptionsfaktoren und ähnliches. Das klingt nicht allzu spezifisch für die Niere. Vielleicht kommen wir später darauf zurück.

### Normalisierung

Zurück zu unserer großen Count-Matrix. Hier nochmal die Ecke links oben.

```{r}
mouse_counts[1:5,1:5]
```

Wie viele Reads haben wir insgesamt in jeder Probe?

Für solche Fragen haben wir zwei Möglichkeiten: Die häufiger verwendete ist, R anzuweisen, die Tabelle als Matrix zu betrachten:

```{r}
mouse_counts %>% column_to_rownames( "gene_id" ) %>% as.matrix() -> count_matrix

count_matrix[ 1:5, 1:5 ]
```

Das hat scheinbar nichts geändert, aber nun stehen uns spezielle Matrix-Funktionen zur Verfügung, z.B. `colSum`, um die Spaltensummen ("column sums") zu berechnen

```{r}
colSums( count_matrix ) %>% head()
```

Wenn wir aber keine neuen Funktionen lernen möchten, können wir auch bei Tidyverse bleiben. Dann müssen wir aber zunächst unsere "breite" Tabelle in eine lange umwandeln, um besser damit arbeiten zu können:

```{r}
mouse_counts %>%
pivot_longer( -gene_id, names_to="sample", values_to="count" ) -> counts_long

head( counts_long )
```

Nun können wir sie Read-Zahl-Summen wie folgt ermitteln:

```{r}
counts_long %>%
group_by( sample ) %>%
summarize( sum(count) ) %>%
head()
```

Wir sehen, dass sich die Gesamtzahl der Reads recht stark von Probe zu Probe variiert. Dies hat aber keine biologische Bedeutung; es spiegelt schlicht wieder, wie gut es gelungen ist, die DNA dieser Probe effizient in die Flowcell zu füllen. Von Bedeutung ist daher nicht die absolute Anzahl der reads, die auf ein Gen gefallen sind, sondern deren *Anteil* an allen Reads der Probe. 

Daher sollten wir jede Read-Zahl durch diese Gesamtzahl teilen. Einfach nur, um bequemere Zahlen zu haben, multiplizieren wir anschließen mit 1 Milliion (`1e6`). Das ergebnis nennen wir "counts per million" ("CPM"):

```{r}
counts_long %>%
group_by( sample ) %>%
mutate( cpm = count / sum(count) * 1e6 ) -> expr_long

head(expr_long)
```

Damit sind nun Werte von verschiedenen Proben miteinander vergleichbar.

Nur als Ergänzung: Wenn man im Matrix-Bild arbeiten möchte, erreicht man dasselbe wie folgt:

```{r}
t( t(count_matrix) / colSums(count_matrix) ) * 1e6 -> cpm_matrix

cpm_matrix[1:5,1:5]
```
Hier haben wir jede Zeile der Maztrix durch ihre Spaltensumme (aus `colSums`) geteilt. Da aber R "matrix / vector" interpretiert als Aufforderung, jede *Zeile* (und nicht: Spalte) der Matrix durch ein Element des vektors zu teilen, müssen wir die Matrix vorher mit `t` transponieren (d.h. Zeilen und Spalten vertauschen) und nachd er Division wieder zurück transponieren. 

Ob man auf diese Weise arbeitet, oder mit Tidyverse, ist Geschmackssache.

### Entwicklung der Leber-Gene: Linienplots

Sehen wir uns nochmal die Liste der Gene an, die in der ersten P0-Leber-Probe am stärsten waren:

```{r}
expr_long %>%
filter( sample == "Liver_P0_1" ) %>%
arrange( desc(count) ) %>%
left_join( gene_names ) %>%
head(10)
```

Wir speichern die IDs dieser 10 Gene zunächst in einem Vektor:

```{r}
expr_long %>%
filter( sample == "Liver_P0_1" ) %>%
arrange( desc(count) ) %>%
left_join( gene_names ) %>%
head(10) %>%
pull( gene_id ) -> top_genes

top_genes
```

Sind diese Gene in den anderen Zeitpunkten auch so hoch? Um dies zu untersuchen, filtern wir die lange Tabelle herunter auf diese Gene, sowie auf alle Leber-Proben:

```{r}
expr_long %>%
filter( gene_id %in% top_genes ) %>%
left_join( sample_table ) %>%
filter( tissue == "Liver" ) %>%
left_join( gene_names ) -> liver_top

head( liver_top )
```

Nun können wir diese Tabelle plotten:

```{r}
liver_top %>%
ggplot( ) +
  geom_point( aes( x=timepoint, y=cpm, col=gene_name ) ) +
  scale_y_log10()
```
Vielleicht wäre es übersichtlicher, wenn wir Linien statt Punkten verwenden würden. Es macht aber keinen Sinn, alle Punkte, die dieselbe Replikat-Ummer haben, mit einer Linie zu verbinden, denn das sind ja dennoch verschiedene Proben. Wir versuchen, über die Replikate zu mitteln:

```{r}
liver_top %>%
group_by( timepoint, gene_name ) %>%
summarise( mean_cpm = mean( cpm ) ) %>%
ggplot( ) +
  geom_line( aes( x=timepoint, y=mean_cpm, group=gene_name, col=gene_name ) ) +
  scale_y_log10()
```

Wir erkennen, dass die meisten der 10 Gene Anfangs schwach exprimiert sind, ab E12.5 schon recht stark und dann langsam weiter zunehmen. Nach der Geburt bleiben sie in etwa konstant. Ein Gen ("Afp") fällt wieder stark ab. "Afp" steht für "alpha-1-Fetoprotein", und der Name zeigt schon, dass es sich um ein Protein handelt, dass wohl nur in Föten vorkommt.

### Log-Expression

Wir haben eben die y-Achse logarithmiert -- und erst so wird der Plot gut lesbar. Das liegt daran, dass Gen-Expressionen über viele Größenordnungen schwanken. 

Nehmen wir dazu nochmal eine Probe heraus, wieder "Liver_P0_1", und plotten ein Histogramm aller cpm-Werte (die wir der entsprechenden Spalte der cpm-Matrix von oben entnehmen):

```{r}
hist( cpm_matrix[,"Liver_P0_1"], 100 )
```
Das war nicht hilfreich.

Wenn wir aber den Logarithmus ziehen, wird es besser:

```{r}
hist( log10( cpm_matrix[,"Liver_P0_1"] ), 100 )
```
Wir sehen, dass die meisten Gene zwischen $10^{0.5}\approx 3$ und $10^2=100$ cpm haben, einige wenige aber über $10^4=10000$.

Klar ist auch, dass die log10(CPM)-Werte, die von -1 bis 5 gehen, besser geeignet scheinen, um damit zu arbeiten, als die CPM-Werte, die zwar theoretisch von 0 bis über 100.000 gehen, meist aber unter 100 bleiben. Deshalb werden wir im folgenden meist mit logarithmierten Werten arbeiten. 

Der Logarithmus von 0 ist minus unendlich, und daher gehen uns alle Nullen in der CPM-Matrix verloren. Um sie zu behalten, zählen wir zu jedem CPM-Wert den Wert 0.03 hinzu. Damit wird 0 zu $\log_{10}(0.03)=-1.5$, und die anderen Werte schieben ein ganz kleines bisschen nach rechts:

```{r}
hist( log10( cpm_matrix[,"Liver_P0_1"] + 0.03 ), 100 )
```

Wie man sieht, haben wir sehr viele Nullen.

Wenn wir nun auf unseren Linien-Plot oben zurück gehen, stellt sich eien Frage. Dort haben wir zwar mit der logarithmierten y-Achse die Expression bereits logarithmisch dargestellt, aber: Sollten wir erst den Logarithmus ziehen und dann über die Replikate mitteln, oder anders herum? Letzteres ist tatsächlich besser, weil sonst ein einzelnes stark exprimiertes Replikat zu viel Einfluss auf den Mittelwert hätte.

Also, hier der Linienplot nochmal, nun mit Mittelung nach dem Logarithmus. 

Wir erstellen erst eine Tabelle mit allen Genen, in der wir den Logarithmus der CPM-Werte nehmen und dann über die Replikate mitteln:

```{r}
expr_long %>%
left_join( sample_table ) %>%
mutate( lcpm = log10( cpm + 0.03 ) ) %>%
group_by( tissue, timepoint, gene_id ) %>%
summarise( mean_lcpm = mean( lcpm ) ) %>%
left_join( gene_names )-> lcpm_means

head( lcpm_means )
```
Nun der Linienplot

```{r}
lcpm_means %>%
filter( tissue == "Liver", gene_name %in% liver_top$gene_name ) %>%
ggplot( ) +
  geom_line( aes( x=timepoint, y=mean_lcpm, group=gene_name, col=gene_name ) ) 
```

### Stark variierende Gene

Unsere Auswahl der 10 Gene war nicht wirklich geeignet, um Gene zu finden, die sich in der Entwicklung stark ändern. Versuchen wir etwas besseres: WIr berechnen für jedes Gen den kleinsten und den größten Wert, den es über alle Zeitpunkte annimmt:

```{r}
lcpm_means %>%
filter( tissue == "Liver" ) %>%
group_by( gene_id, gene_name ) %>%
summarise( 
  min_lcpm = min( mean_lcpm ),
  max_lcpm = max( mean_lcpm ) ) %>%
arrange( desc( max_lcpm - min_lcpm ) ) %>%
head()
```

Beachten Sie: Wir haben im vorigen Schritt über die Replikate gemittelt. Nun haben wir alle diese Mittelwerte nach Genen gruppiert, d.h., die Werte vom selben Gen, aber von verschiedenen Zeitpunkten zusammen genommen, und hier die Minima und Maxima ermittelt. Dann haben wir nach der Differenz zwischen Maximum und Minimum (also der Spanne der Werte) absteigend sortiert.

Wir schreiben uns die 50 Gene mit der größten Spanne heraus:

```{r}
lcpm_means %>%
filter( tissue == "Liver" ) %>%
group_by( gene_id, gene_name ) %>%
summarise( 
  min_lcpm = min( mean_lcpm ),
  max_lcpm = max( mean_lcpm ) ) %>%
arrange( desc( max_lcpm - min_lcpm ) ) %>%

head( 50 ) %>%
pull( gene_id ) -> top_50_by_span

head( top_50_by_span )
```

Der Linienplot für diese Gene sieht nun etwas unübersichtlich aus:

```{r}
lcpm_means %>%
filter( tissue == "Liver", gene_id %in% top_50_by_span ) %>%
ggplot( ) +
  geom_line( aes( x=timepoint, y=mean_lcpm, group=gene_name, col=gene_name ) ) 
```
Wie können wir dies besser darstellen?

### Heatmaps

Hier ist eine andere Darstellung derselben Daten, nämlich mit `geom_tile` statt `geom_line`. Wir tauschen y-Achse gegen Farbe: die gene unterscheiden wir nun, indem wir sie untereinander anordnen, und die Expressionsstärke stellen wir durch die Farbe dar.

```{r fig.height=11, fig.width=8}
lcpm_means %>%
filter( tissue == "Liver", gene_id %in% top_50_by_span ) %>%
ggplot( ) +
  geom_tile( aes( x=timepoint, y=gene_name, fill=mean_lcpm ) ) +
  scale_fill_viridis_c( option="magma" )
```

Etwas unübersichtlich ist es immer noch. Es hilt, die Gene anders zu sortieren. 

Es gibt spezielel Funktionen um in solchen Heatmaps die Zeilen doer Spalten so zu sortieren, dass man besser sehen kann, was sich ähnelt. Wir werden das Paket `pheatmap` (pretty heatmaps) benutzen:

```{r}
library( pheatmap )
```

Während wir ggplot mit geom_tile eben die daten in einer langen Liste übergeben haben, möchte pheatmap sie als Matrix. Die Zeitpunkte sollten also nicht untereinander, sondern nebeneinander stehen. Wir benutzen `pivot_wider` hierzu, und wandeln dann noch alles zu einer Matrix um, die die Gen-Namen als Zeilen-Namen hat:

```{r}
lcpm_means %>%
filter( tissue == "Liver", gene_id %in% top_50_by_span ) %>%
pivot_wider( id_cols = "gene_name", names_from = "timepoint", values_from = "mean_lcpm" ) %>%
column_to_rownames( "gene_name" ) %>%
as.matrix() -> liver_top50_varying_matrix

liver_top50_varying_matrix[ 1:5, 1:7 ]
```

Nun rufen wir die Funktion `pheatmap` aus dem paket `pheatmap` auf:


```{r fig.height=11, fig.width=8}
pheatmap( 
  liver_top50_varying_matrix,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  clustering_distance_rows = "correlation",
  color = viridis::magma(100) )
```

Die Funktion `pheatmap` hat Dutzende Argumente, um das Aussehen der Heatmap anzupassen. Wir haben hier nur einige verwendet.

Das erste Argument ist die Matrix mit den Daten. Mit `cluster_cols=FALSE` und `cluster_rows=TRUE` stellen wir klar, dass die Reihenfolge der Spalten unverändert aus der Matrix übernommen werden sollen, während die Spalten durch sog. Clustering umgeordnet werden sollen. Zuletzt wählen wir mit `color = viridis::magma(100)` die Farbskala aus: Wir möchten die "Magma"-Palette aus dem "viridis"-Paket in 100 Abstufungungen verwenden. Mehr zu den Viridis-Paletten fidnen Sie [hier](https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html).

Beim Clustering wird zunächst bestimmt, wie "ähnlich" die Zeilen zueinander sind: Zwei Zeilen gelten als "nah beieinander", wenn sie einen hohen Korrelationskoeffizienten haben. (Zum Korrelationskoeffizienten später mehr.) Dass zur Beurteilung der Ähnlichkeit der Korrelationskoeffizient verwendet werden soll, wird durch `clustering_distance_rows = "correlation"` festgelegt. 

Hier wird sog. hierarchisches (oder: agglomeratives) Clustering verwendet, um Gruppen ähnlicher Zielen nebeneinander zu setzen. Näheres heirzu findet sich [hier](https://de.wikipedia.org/wiki/Hierarchische_Clusteranalyse) auf Wikipedia.

#### Beurteilung

Mit der verbesserten Heatmap können wir jetzt einiges ablesen: Die Gene, die nur im Fötus verwendet werden, haben alle mit fötalem Hämoglobin zu tun. (Erinnern Sie sich, die Blutbildung im erwachsenen Säugetier zwar nur im Knochenmark erfolgt, im Embryo aber auch in Leber und Milz, und im Fötus ein anderes Hämoglobin als im ausgewachsens Tier verwendet wird. Siehe auch [hier](https://en.wikipedia.org/wiki/Fetal_hemoglobin).)

Wir können auch Gruppen von Genen erkennen, die erst im reiferen Fötus aktiv werden, und einige, die sogar erst im ausgewachsenen Tier aktiv werden. Sie werdn sicher erkennen, dass wir viel über die Leber lernen können, wenn wir solche Plots im Detail studieren.

### Reproduzierbarkeit

Erstaunlich ist das Gen Sult2a1, das scheinbar nur zu einem Zeitpunkt, nämlich 14 Tage nach Geburt (P14) hoch exprimiert ist, und dann wieder abnimmt. Stimmt das, oder ist hier vielleicht ein einzelnes der 4 Replikate ein Ausreißer und hat den Mittwelwert für diesen Zeitpunkt hoch gezogen?

Wir prüfen das, indem wir auf die ursprünglichen, noch nicht über Replikate gemittelten Daten zurück gehen und das Gen plotten:

```{r}
expr_long %>%
left_join( sample_table ) %>%
filter( tissue == "Liver" ) %>%
left_join( gene_names ) %>%
filter( gene_name == "Sult2a1" ) %>%
ggplot + 
  geom_point( aes( x=timepoint, y=cpm ) ) +
  scale_y_log10()
```

Wir erkennen: Das das Gen langsam ansteigt, und zum Zeitpunkt P14 maximal ist, da sind sich alle Replikate einig, auch zu dem Abfall kurz vor Geburt. Der Effekt im höheren Alter ist aber zwischen Replikaten sehr variabel.

Wir möchten aber auch nicht zu viele Plots ansehen müssen. Daher gibt es statistsiche Techniken, um automatisiert zu beurteilen, für welche Gene Effekte als signifikant gesehen werden dürfen, da alle Repplikate ein konsitentes Bild zeigen, und wo die Daten unklarer sind. Dies wird als "differential expression calling" oder als Ermittlung von "differentially expressed genes" (DGE) bezeichnet. Bei einem einfachen Vergleich zwischen zwei Gruppen handelt es sich dabei um Weiterentwicklungen aus dem t-Test heraus. Für kompliziertere Situationen, wie die Zeitreihe hier, ist dies aber etwas komplizierter. Eventuell werden wir in einer späteren Vorlesung auf Details eingehen können.

### Proben-Korrelation

Die Frage, ob die Replikate gut miteinander übereinstimmen, können wir auch beurteilen, indem wir sie gegeneinander plotten. Wir verwenden hierzu unsere Expressionsdaten in Matrixform (siehe oben), logarithmiert:

```{r}
log10( cpm_matrix + .03 )[ 1:5, 1:5 ]
```

Nun plotten wir zwei Spalten gegeneinander, z.B. die ersten beiden Leber-Proben zu P0:

```{r}
plot( 
  log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
  log10( cpm_matrix + .03 )[ , "Liver_P0_2" ],
  cex=.2, asp=1, col = alpha( "black", .3 ) )
abline( 0, 1, col="orange")
```
(Zu meiner Bequenlichkeit habe ich hier die `plot`-Funktion aus Basis-R, nicht die aus ggplot, verwendet.)

Zum Vergleich, zwei Proben von verschiedenen Zeitpunkten:

```{r}
plot( 
  log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
  log10( cpm_matrix + .03 )[ , "Liver_e15.5_1" ],
  cex=.2, asp=1, col = alpha( "black", .3 ) )
abline( 0, 1, col="orange")
```

Wir erkennen stärkere Abweichung. 

Wir können dies auch durch den Korrelationskoeffizienten quantifizieren, der im zweiten Fall niedriger ist.

```{r}
cor( 
  log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
  log10( cpm_matrix + .03 )[ , "Liver_P0_2" ] )

cor( 
  log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
  log10( cpm_matrix + .03 )[ , "Liver_e15.5_1" ] )
```

Man kann der Funktion `cor` auch statt zweier Spalten die ganze Matrix übergeben, und erhält dann eine sog. Korrelationsmatrix mit der Korrelation aller Spalten zu allen Spalten. Diese matrxi können wir dann mit `pheatmap` plotten:

```{r fig.height=10, fig.width=9}
pheatmap(
  cor( log10( cpm_matrix + .03 ) ),
  cluster_rows=FALSE, cluster_cols=FALSE, 
  annotation_col = sample_table %>% select( - replicate ) %>% column_to_rownames( "sample" ),
  labels_row = "", labels_col = "" )
```

Hier dasselbe, nur für die Leber-Proben:

```{r fig.height=10, fig.width=9}
pheatmap(
  cor( log10( cpm_matrix[ , sample_table$tissue=="Liver" ] + .03 ) ),
  cluster_rows=FALSE, cluster_cols=FALSE, 
  annotation_col = sample_table %>% select( - replicate ) %>% column_to_rownames( "sample" ) )
```
Hieraus können wir nun ablesen, zu welchen Zeitpunkten sich das Leber-Transkriptom besonder stark ändert.
 
Aber wir erkennen auch Probleme mit der Probenqualität. Einige Proben unterscheiden sich von ihren replikaten deutlich stärker als anders, z.B. die allererste Probe.
 
 
### Gewebespezifische Gene

Weiter oben haben wir gefragt, welche Gene in der Niere besodners stark exprimiert sind, konnten aber mit den gefunden Genen nicht viel anfangen. Daher stellen wir die Frage nun etwas anders:

Welche Gene sind besonders spezifisch für die Niere? Wir werden diese Frage für ausgewachsene Mäuse beantworten, indem wir uns auf die Proben von Zeitpunkt P28 beschränken.

Zunächst nehmen wir die vier P28-Proben von der Niere und bestimmen für jedes Gen den jeweils kleinsten der vier CPM-Werte:

```{r}
expr_long %>%
left_join( sample_table ) %>%
filter( timepoint == "P28", tissue == "Kidney" ) %>%
group_by( gene_id ) %>%
summarise( min_cpm_in_kidney = min( cpm ) ) -> min_kidney_P28

min_kidney_P28
```

Nun bestimmen wir für die P28-Proben von allen Geweben *außer* der Niere den jeweils größten Wert:

```{r}
expr_long %>%
left_join( sample_table ) %>%
filter( timepoint == "P28", tissue != "Kidney" ) %>%
group_by( gene_id ) %>%
summarise( max_cpm_outside_kidney = max( cpm ) ) -> max_non_kidney_P28

max_non_kidney_P28
```

Wenn wir diese beiden Tabellen nun nebeneinander setzen, können wir die Gene bestimmen, die in der Niere viel höher sind als in all den anderen Proben:

```{r}
inner_join( min_kidney_P28, max_non_kidney_P28 ) %>%
arrange( desc( min_cpm_in_kidney / ( max_cpm_outside_kidney + 0.03 ) ) ) %>%
left_join( gene_names )
```
```{r}
inner_join(
  
  lcpm_means %>%
  filter( timepoint == "P28", tissue == "Kidney" ) %>%
  rename( kidney_lcpm = mean_lcpm ),

  lcpm_means %>%
  filter( timepoint == "P28", tissue != "Kidney" ) %>%
  group_by( gene_name ) %>%
  summarise( non_kidney_max_lcpm = max( mean_lcpm ) ) ) %>%

mutate( kidney_lead = kidney_lcpm - non_kidney_max_lcpm ) %>%
arrange( desc( kidney_lead ) )
```