Hier verwenden wir Daten aus folgendem Paper:
M. Cardoso-Moreia et al.: Gene expression across mammalian organ development. Nature 571:505 (2019).
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 für Details); die Matrix, die ich so erhalten habe, ist hier hier verfügbar).
Wir laden diese Matrix
suppressPackageStartupMessages( library(tidyverse) )
read_tsv( "Downloads/evodevo_mouse_counts.tsv.gz" ) -> mouse_countsRows: 36141 Columns: 317
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): gene_id
dbl (316): Brain_e10.5_1, Brain_e10.5_2, Brain_e10.5_3, Brain_e11.5_1, Brain...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Hier sind die ersten Spalten und Zeilen der Matrix:
mouse_counts[ 1:5, 1:5 ]# A tibble: 5 × 5
gene_id Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
<chr> <dbl> <dbl> <dbl> <dbl>
1 ENSMUSG00000000001 9532 7451 7631 8833
2 ENSMUSG00000000003 0 0 0 0
3 ENSMUSG00000000028 2374 2299 2112 1724
4 ENSMUSG00000000037 322 282 263 328
5 ENSMUSG00000000049 0 0 0 0Die 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.
tibble( sample = colnames(mouse_counts)[-1] ) %>% head()# A tibble: 6 × 1
sample
<chr>
1 Brain_e10.5_1
2 Brain_e10.5_2
3 Brain_e10.5_3
4 Brain_e11.5_1
5 Brain_e11.5_2
6 Brain_e11.5_3Mit dem Tidyverse-Verb separate können wir die Inhalte der sample-Spalte entlang der Unterstriche aufteilen und in drei neue Spalten schreiben:
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 )# A tibble: 6 × 4
sample tissue timepoint replicate
<chr> <chr> <fct> <chr>
1 Brain_e10.5_1 Brain e10.5 1
2 Brain_e10.5_2 Brain e10.5 2
3 Brain_e10.5_3 Brain e10.5 3
4 Brain_e11.5_1 Brain e11.5 1
5 Brain_e11.5_2 Brain e11.5 2
6 Brain_e11.5_3 Brain e11.5 3 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:
sample_table %>% pull(timepoint) %>% levels() [1] "e10.5" "e11.5" "e12.5" "e13.5" "e14.5" "e15.5" "e16.5" "e17.5" "e18.5"
[10] "P0" "P3" "P14" "P28" "P63" 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:
sample_table %>% pull( tissue ) %>% unique()[1] "Brain" "Cerebellum" "Heart" "Kidney" "Liver"
[6] "Ovary" "Testis" (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:
sample_table %>%
group_by( tissue, timepoint ) %>%
summarise( n=n() ) %>%
pivot_wider( names_from = "tissue", values_from = "n", values_fill=0 )`summarise()` has grouped output by 'tissue'. You can override using the
`.groups` argument.# A tibble: 14 × 8
timepoint Brain Cerebellum Heart Kidney Liver Ovary Testis
<fct> <int> <int> <int> <int> <int> <int> <int>
1 e10.5 3 0 4 0 4 2 2
2 e11.5 4 0 4 4 4 2 2
3 e12.5 4 0 4 4 4 2 2
4 e13.5 4 4 4 4 5 2 2
5 e14.5 4 4 3 3 4 2 2
6 e15.5 4 4 4 4 4 2 2
7 e16.5 4 4 4 4 4 2 2
8 e17.5 4 4 4 4 4 2 2
9 e18.5 4 4 4 4 4 2 1
10 P0 4 4 4 4 4 2 2
11 P3 4 4 4 4 4 2 2
12 P14 4 3 4 4 4 2 2
13 P28 4 4 4 4 4 2 2
14 P63 4 4 4 4 4 2 2Erkennen 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?
Nun sollten wir einen ersten Blick auf unsere Count-Matrix werfen:
mouse_counts[1:5,1:5]# A tibble: 5 × 5
gene_id Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
<chr> <dbl> <dbl> <dbl> <dbl>
1 ENSMUSG00000000001 9532 7451 7631 8833
2 ENSMUSG00000000003 0 0 0 0
3 ENSMUSG00000000028 2374 2299 2112 1724
4 ENSMUSG00000000037 322 282 263 328
5 ENSMUSG00000000049 0 0 0 0Welche 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”)
mouse_counts %>%
select( gene_id, Liver_P0_1 ) %>%
arrange( desc( Liver_P0_1 ) )# A tibble: 36,141 × 2
gene_id Liver_P0_1
<chr> <dbl>
1 ENSMUSG00000029368 1608170
2 ENSMUSG00000022868 496738
3 ENSMUSG00000024164 247633
4 ENSMUSG00000054932 238008
5 ENSMUSG00000032083 195283
6 ENSMUSG00000002985 139282
7 ENSMUSG00000028001 134252
8 ENSMUSG00000030359 127224
9 ENSMUSG00000033831 115231
10 ENSMUSG00000064373 110436
# ℹ 36,131 more rowsMit 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 erstellt haben. Wir benennen aber die Spalten um, um sie etwas kürzer zu machen:
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_namesRows: 56305 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): Gene stable ID, Gene name, Gene description
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head( gene_names )# A tibble: 6 × 3
gene_id gene_name gene_descr
<chr> <chr> <chr>
1 ENSMUSG00000064372 mt-Tp mitochondrially encoded tRNA proline [Source:MGI…
2 ENSMUSG00000064371 mt-Tt mitochondrially encoded tRNA threonine [Source:M…
3 ENSMUSG00000064370 mt-Cytb mitochondrially encoded cytochrome b [Source:MGI…
4 ENSMUSG00000064369 mt-Te mitochondrially encoded tRNA glutamic acid [Sour…
5 ENSMUSG00000064368 mt-Nd6 mitochondrially encoded NADH dehydrogenase 6 [So…
6 ENSMUSG00000064367 mt-Nd5 mitochondrially encoded NADH dehydrogenase 5 [So…Wir fügen diese Informationen unserem Ergebnis an, indem wir den Code von eben durch ein left_join ergänzen:
mouse_counts %>%
select( gene_id, Liver_P0_1 ) %>%
arrange( desc( Liver_P0_1 ) ) %>%
left_join( gene_names ) %>%
head( 10 )Joining with `by = join_by(gene_id)`# A tibble: 10 × 4
gene_id Liver_P0_1 gene_name gene_descr
<chr> <dbl> <chr> <chr>
1 ENSMUSG00000029368 1608170 Alb albumin [Source:MGI Symbol;Acc:MGI:8…
2 ENSMUSG00000022868 496738 Ahsg alpha-2-HS-glycoprotein [Source:MGI …
3 ENSMUSG00000024164 247633 C3 complement component 3 [Source:MGI S…
4 ENSMUSG00000054932 238008 Afp alpha fetoprotein [Source:MGI Symbol…
5 ENSMUSG00000032083 195283 Apoa1 apolipoprotein A-I [Source:MGI Symbo…
6 ENSMUSG00000002985 139282 Apoe apolipoprotein E [Source:MGI Symbol;…
7 ENSMUSG00000028001 134252 Fga fibrinogen alpha chain [Source:MGI S…
8 ENSMUSG00000030359 127224 Pzp PZP, alpha-2-macroglobulin like [Sou…
9 ENSMUSG00000033831 115231 Fgb fibrinogen beta chain [Source:MGI Sy…
10 ENSMUSG00000064373 110436 Selenop selenoprotein P [Source:MGI Symbol;A…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:
mouse_counts %>%
select( gene_id, Kidney_P0_1 ) %>%
arrange( desc( Kidney_P0_1 ) ) %>%
left_join( gene_names ) %>%
head( 10 )Joining with `by = join_by(gene_id)`# A tibble: 10 × 4
gene_id Kidney_P0_1 gene_name gene_descr
<chr> <dbl> <chr> <chr>
1 ENSMUSG00000037742 94461 Eef1a1 eukaryotic translation elongation f…
2 ENSMUSG00000033161 70392 Atp1a1 ATPase, Na+/K+ transporting, alpha …
3 ENSMUSG00000031502 61047 Col4a1 collagen, type IV, alpha 1 [Source:…
4 ENSMUSG00000034994 49973 Eef2 eukaryotic translation elongation f…
5 ENSMUSG00000023944 48543 Hsp90ab1 heat shock protein 90 alpha (cytoso…
6 ENSMUSG00000032216 47811 Nedd4 neural precursor cell expressed, de…
7 ENSMUSG00000025393 41872 Atp5b ATP synthase, H+ transporting mitoc…
8 ENSMUSG00000004980 39683 Hnrnpa2b1 heterogeneous nuclear ribonucleopro…
9 ENSMUSG00000022816 38032 Fstl1 follistatin-like 1 [Source:MGI Symb…
10 ENSMUSG00000029580 37909 Actb actin, beta [Source:MGI Symbol;Acc:…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.
Zurück zu unserer großen Count-Matrix. Hier nochmal die Ecke links oben.
mouse_counts[1:5,1:5]# A tibble: 5 × 5
gene_id Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
<chr> <dbl> <dbl> <dbl> <dbl>
1 ENSMUSG00000000001 9532 7451 7631 8833
2 ENSMUSG00000000003 0 0 0 0
3 ENSMUSG00000000028 2374 2299 2112 1724
4 ENSMUSG00000000037 322 282 263 328
5 ENSMUSG00000000049 0 0 0 0Wie 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:
mouse_counts %>% column_to_rownames( "gene_id" ) %>% as.matrix() -> count_matrix
count_matrix[ 1:5, 1:5 ] Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
ENSMUSG00000000001 9532 7451 7631 8833
ENSMUSG00000000003 0 0 0 0
ENSMUSG00000000028 2374 2299 2112 1724
ENSMUSG00000000037 322 282 263 328
ENSMUSG00000000049 0 0 0 0
Brain_e11.5_2
ENSMUSG00000000001 11079
ENSMUSG00000000003 0
ENSMUSG00000000028 2064
ENSMUSG00000000037 371
ENSMUSG00000000049 0Das 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
colSums( count_matrix ) %>% head()Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1 Brain_e11.5_2
28702741 28346753 28061985 25885460 29435786
Brain_e11.5_3
30547013 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:
mouse_counts %>%
pivot_longer( -gene_id, names_to="sample", values_to="count" ) -> counts_long
head( counts_long )# A tibble: 6 × 3
gene_id sample count
<chr> <chr> <dbl>
1 ENSMUSG00000000001 Brain_e10.5_1 9532
2 ENSMUSG00000000001 Brain_e10.5_2 7451
3 ENSMUSG00000000001 Brain_e10.5_3 7631
4 ENSMUSG00000000001 Brain_e11.5_1 8833
5 ENSMUSG00000000001 Brain_e11.5_2 11079
6 ENSMUSG00000000001 Brain_e11.5_3 11377Nun können wir sie Read-Zahl-Summen wie folgt ermitteln:
counts_long %>%
group_by( sample ) %>%
summarize( sum(count) ) %>%
head()# A tibble: 6 × 2
sample `sum(count)`
<chr> <dbl>
1 Brain_P0_1 44411177
2 Brain_P0_2 16292683
3 Brain_P0_3 13873908
4 Brain_P0_4 16032449
5 Brain_P14_1 14312695
6 Brain_P14_2 16828137Wir 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”):
counts_long %>%
group_by( sample ) %>%
mutate( cpm = count / sum(count) * 1e6 ) -> expr_long
head(expr_long)# A tibble: 6 × 4
# Groups: sample [6]
gene_id sample count cpm
<chr> <chr> <dbl> <dbl>
1 ENSMUSG00000000001 Brain_e10.5_1 9532 332.
2 ENSMUSG00000000001 Brain_e10.5_2 7451 263.
3 ENSMUSG00000000001 Brain_e10.5_3 7631 272.
4 ENSMUSG00000000001 Brain_e11.5_1 8833 341.
5 ENSMUSG00000000001 Brain_e11.5_2 11079 376.
6 ENSMUSG00000000001 Brain_e11.5_3 11377 372.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:
t( t(count_matrix) / colSums(count_matrix) ) * 1e6 -> cpm_matrix
cpm_matrix[1:5,1:5] Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
ENSMUSG00000000001 332.09372 262.851975 271.93372 341.23404
ENSMUSG00000000003 0.00000 0.000000 0.00000 0.00000
ENSMUSG00000000028 82.70987 81.102763 75.26196 66.60110
ENSMUSG00000000037 11.21844 9.948229 9.37211 12.67121
ENSMUSG00000000049 0.00000 0.000000 0.00000 0.00000
Brain_e11.5_2
ENSMUSG00000000001 376.37860
ENSMUSG00000000003 0.00000
ENSMUSG00000000028 70.11873
ENSMUSG00000000037 12.60371
ENSMUSG00000000049 0.00000Hier 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.
Sehen wir uns nochmal die Liste der Gene an, die in der ersten P0-Leber-Probe am stärsten waren:
expr_long %>%
filter( sample == "Liver_P0_1" ) %>%
arrange( desc(count) ) %>%
left_join( gene_names ) %>%
head(10)Joining with `by = join_by(gene_id)`# A tibble: 10 × 6
# Groups: sample [1]
gene_id sample count cpm gene_name gene_descr
<chr> <chr> <dbl> <dbl> <chr> <chr>
1 ENSMUSG00000029368 Liver_P0_1 1608170 118978. Alb albumin [Source:MGI …
2 ENSMUSG00000022868 Liver_P0_1 496738 36750. Ahsg alpha-2-HS-glycoprot…
3 ENSMUSG00000024164 Liver_P0_1 247633 18321. C3 complement component…
4 ENSMUSG00000054932 Liver_P0_1 238008 17609. Afp alpha fetoprotein [S…
5 ENSMUSG00000032083 Liver_P0_1 195283 14448. Apoa1 apolipoprotein A-I […
6 ENSMUSG00000002985 Liver_P0_1 139282 10305. Apoe apolipoprotein E [So…
7 ENSMUSG00000028001 Liver_P0_1 134252 9932. Fga fibrinogen alpha cha…
8 ENSMUSG00000030359 Liver_P0_1 127224 9412. Pzp PZP, alpha-2-macrogl…
9 ENSMUSG00000033831 Liver_P0_1 115231 8525. Fgb fibrinogen beta chai…
10 ENSMUSG00000064373 Liver_P0_1 110436 8170. Selenop selenoprotein P [Sou…Wir speichern die IDs dieser 10 Gene zunächst in einem Vektor:
expr_long %>%
filter( sample == "Liver_P0_1" ) %>%
arrange( desc(count) ) %>%
left_join( gene_names ) %>%
head(10) %>%
pull( gene_id ) -> top_genesJoining with `by = join_by(gene_id)`top_genes [1] "ENSMUSG00000029368" "ENSMUSG00000022868" "ENSMUSG00000024164"
[4] "ENSMUSG00000054932" "ENSMUSG00000032083" "ENSMUSG00000002985"
[7] "ENSMUSG00000028001" "ENSMUSG00000030359" "ENSMUSG00000033831"
[10] "ENSMUSG00000064373"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:
expr_long %>%
filter( gene_id %in% top_genes ) %>%
left_join( sample_table ) %>%
filter( tissue == "Liver" ) %>%
left_join( gene_names ) -> liver_topJoining with `by = join_by(sample)`
Joining with `by = join_by(gene_id)`head( liver_top )# A tibble: 6 × 9
# Groups: sample [6]
gene_id sample count cpm tissue timepoint replicate gene_name gene_descr
<chr> <chr> <dbl> <dbl> <chr> <fct> <chr> <chr> <chr>
1 ENSMUSG000… Liver… 1426 61.4 Liver e10.5 1 Apoe apolipopr…
2 ENSMUSG000… Liver… 11994 456. Liver e10.5 2 Apoe apolipopr…
3 ENSMUSG000… Liver… 13957 542. Liver e10.5 3 Apoe apolipopr…
4 ENSMUSG000… Liver… 5813 308. Liver e10.5 4 Apoe apolipopr…
5 ENSMUSG000… Liver… 5420 221. Liver e11.5 2 Apoe apolipopr…
6 ENSMUSG000… Liver… 8061 217. Liver e11.5 3 Apoe apolipopr…Nun können wir diese Tabelle plotten:
liver_top %>%
ggplot( ) +
geom_point( aes( x=timepoint, y=cpm, col=gene_name ) ) +
scale_y_log10()Warning: Transformation introduced infinite values in continuous y-axis
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:
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()`summarise()` has grouped output by 'timepoint'. You can override using the
`.groups` argument.
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.
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):
hist( cpm_matrix[,"Liver_P0_1"], 100 )
Das war nicht hilfreich.
Wenn wir aber den Logarithmus ziehen, wird es besser:
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:
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:
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_meansJoining with `by = join_by(sample)`
`summarise()` has grouped output by 'tissue', 'timepoint'. You can override
using the `.groups` argument.
Joining with `by = join_by(gene_id)`head( lcpm_means )# A tibble: 6 × 6
# Groups: tissue, timepoint [1]
tissue timepoint gene_id mean_lcpm gene_name gene_descr
<chr> <fct> <chr> <dbl> <chr> <chr>
1 Brain e10.5 ENSMUSG00000000001 2.46 Gnai3 guanine nucleotide bi…
2 Brain e10.5 ENSMUSG00000000003 -1.52 Pbsn probasin [Source:MGI …
3 Brain e10.5 ENSMUSG00000000028 1.90 Cdc45 cell division cycle 4…
4 Brain e10.5 ENSMUSG00000000037 1.01 Scml2 Scm polycomb group pr…
5 Brain e10.5 ENSMUSG00000000049 -1.52 Apoh apolipoprotein H [Sou…
6 Brain e10.5 ENSMUSG00000000056 1.98 Narf nuclear prelamin A re…Nun der Linienplot
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 ) ) 
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:
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()`summarise()` has grouped output by 'gene_id'. You can override using the
`.groups` argument.# A tibble: 6 × 4
# Groups: gene_id [6]
gene_id gene_name min_lcpm max_lcpm
<chr> <chr> <dbl> <dbl>
1 ENSMUSG00000058207 Serpina3k -1.52 4.29
2 ENSMUSG00000052187 Hbb-y -1.52 3.83
3 ENSMUSG00000056035 Cyp3a11 -1.52 3.79
4 ENSMUSG00000025479 Cyp2e1 -1.43 3.79
5 ENSMUSG00000066154 Mup3 -1.52 3.65
6 ENSMUSG00000052217 Hbb-bh1 -1.52 3.61Beachten 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:
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`summarise()` has grouped output by 'gene_id'. You can override using the
`.groups` argument.head( top_50_by_span )[1] "ENSMUSG00000058207" "ENSMUSG00000052187" "ENSMUSG00000056035"
[4] "ENSMUSG00000025479" "ENSMUSG00000066154" "ENSMUSG00000052217"Der Linienplot für diese Gene sieht nun etwas unübersichtlich aus:
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?
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.
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:
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:
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 ] e10.5 e11.5 e12.5 e13.5 e14.5 e15.5
Tat -1.095259 -1.276733 -0.5821924 0.08844903 0.0181510 0.2505338
Cyp2c29 -1.433955 -1.522879 -1.5228787 -1.52287875 -1.4192102 -1.5228787
Cyp2a5 -1.522879 -1.522879 -1.2586920 -0.49534673 -0.4743501 0.1250718
Slc10a1 -1.522879 -1.429619 -0.6561367 0.42552969 0.4621024 1.2915968
Akr1c6 -1.150737 -1.522879 -0.3346262 0.59366773 0.9175212 0.9349686
e16.5
Tat 0.5080352
Cyp2c29 -1.3783121
Cyp2a5 0.6813330
Slc10a1 1.7939137
Akr1c6 1.4435632Nun rufen wir die Funktion pheatmap aus dem paket pheatmap auf:
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.
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 auf Wikipedia.
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.)
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.
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:
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()Joining with `by = join_by(sample)`
Joining with `by = join_by(gene_id)`Warning: Transformation introduced infinite values in continuous y-axis
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.
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:
log10( cpm_matrix + .03 )[ 1:5, 1:5 ] Brain_e10.5_1 Brain_e10.5_2 Brain_e10.5_3 Brain_e11.5_1
ENSMUSG00000000001 2.521300 2.4197608 2.4345110 2.533091
ENSMUSG00000000003 -1.522879 -1.5228787 -1.5228787 -1.522879
ENSMUSG00000000028 1.917715 1.9091963 1.8767486 1.823677
ENSMUSG00000000037 1.051092 0.9990535 0.9732253 1.103845
ENSMUSG00000000049 -1.522879 -1.5228787 -1.5228787 -1.522879
Brain_e11.5_2
ENSMUSG00000000001 2.575660
ENSMUSG00000000003 -1.522879
ENSMUSG00000000028 1.846020
ENSMUSG00000000037 1.101531
ENSMUSG00000000049 -1.522879Nun plotten wir zwei Spalten gegeneinander, z.B. die ersten beiden Leber-Proben zu P0:
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:
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.
cor(
log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
log10( cpm_matrix + .03 )[ , "Liver_P0_2" ] )[1] 0.9856564cor(
log10( cpm_matrix + .03 )[ , "Liver_P0_1" ],
log10( cpm_matrix + .03 )[ , "Liver_e15.5_1" ] )[1] 0.9663852Man 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:
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:
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.
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:
expr_long %>%
left_join( sample_table ) %>%
filter( timepoint == "P28", tissue == "Kidney" ) %>%
group_by( gene_id ) %>%
summarise( min_cpm_in_kidney = min( cpm ) ) -> min_kidney_P28Joining with `by = join_by(sample)`min_kidney_P28# A tibble: 36,141 × 2
gene_id min_cpm_in_kidney
<chr> <dbl>
1 ENSMUSG00000000001 170.
2 ENSMUSG00000000003 0
3 ENSMUSG00000000028 3.22
4 ENSMUSG00000000037 0.846
5 ENSMUSG00000000049 0
6 ENSMUSG00000000056 97.7
7 ENSMUSG00000000058 66.1
8 ENSMUSG00000000078 37.4
9 ENSMUSG00000000085 25.9
10 ENSMUSG00000000088 209.
# ℹ 36,131 more rowsNun bestimmen wir für die P28-Proben von allen Geweben außer der Niere den jeweils größten Wert:
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_P28Joining with `by = join_by(sample)`max_non_kidney_P28# A tibble: 36,141 × 2
gene_id max_cpm_outside_kidney
<chr> <dbl>
1 ENSMUSG00000000001 223.
2 ENSMUSG00000000003 0
3 ENSMUSG00000000028 42.2
4 ENSMUSG00000000037 33.4
5 ENSMUSG00000000049 2535.
6 ENSMUSG00000000056 235.
7 ENSMUSG00000000058 138.
8 ENSMUSG00000000078 124.
9 ENSMUSG00000000085 139.
10 ENSMUSG00000000088 680.
# ℹ 36,131 more rowsWenn 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:
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 )Joining with `by = join_by(gene_id)`
Joining with `by = join_by(gene_id)`# A tibble: 36,141 × 5
gene_id min_cpm_in_kidney max_cpm_outside_kidney gene_name gene_descr
<chr> <dbl> <dbl> <chr> <chr>
1 ENSMUSG0000002… 4159. 2.69 Slc34a1 solute ca…
2 ENSMUSG0000004… 564. 0.347 Slc5a12 solute ca…
3 ENSMUSG0000004… 279. 0.174 Lrrc19 leucine r…
4 ENSMUSG0000002… 225. 0.174 Slc22a19 solute ca…
5 ENSMUSG0000002… 1371. 1.48 Mep1a meprin 1 …
6 ENSMUSG0000003… 3417. 3.81 Umod uromoduli…
7 ENSMUSG0000000… 286. 0.336 Dnase1 deoxyribo…
8 ENSMUSG0000004… 371. 0.449 Tmem174 transmemb…
9 ENSMUSG0000002… 1264. 1.78 Slc12a1 solute ca…
10 ENSMUSG0000005… 115. 0.139 Clrn3 clarin 3 …
# ℹ 36,131 more rowsinner_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 ) )Joining with `by = join_by(gene_name)`# A tibble: 36,141 × 8
# Groups: tissue, timepoint [1]
tissue timepoint gene_id kidney_lcpm gene_name gene_descr non_kidney_max_lcpm
<chr> <fct> <chr> <dbl> <chr> <chr> <dbl>
1 Kidney P28 ENSMUS… 4.15 Kap kidney an… 0.727
2 Kidney P28 ENSMUS… 2.48 Lrrc19 leucine r… -0.911
3 Kidney P28 ENSMUS… 2.56 Dnase1 deoxyribo… -0.768
4 Kidney P28 ENSMUS… 2.81 Slc5a12 solute ca… -0.509
5 Kidney P28 ENSMUS… 2.50 Slc22a19 solute ca… -0.811
6 Kidney P28 ENSMUS… 3.66 Slc34a1 solute ca… 0.353
7 Kidney P28 ENSMUS… 3.28 Slc12a1 solute ca… 0.0352
8 Kidney P28 ENSMUS… 2.09 Clrn3 clarin 3 … -1.15
9 Kidney P28 ENSMUS… 2.40 Cyp2j11 cytochrom… -0.843
10 Kidney P28 ENSMUS… 3.41 Napsa napsin A … 0.173
# ℹ 36,131 more rows
# ℹ 1 more variable: kidney_lead <dbl>