Hausaufgabe RNA-Seq

In dieser Hausaufgabe arbeiten wir mit den RNA-Seq-Daten zur Maus aus dem in der Vorlesungen besprochenen Paper von Cardoso-Moreia et al. (2019).

Hier finden Sie eine Datei mit der Count-Matrix, die Sie mit read_tsv laden können.

library( tidyverse ) 

read_tsv( "data_on_git/evodevo_mouse_counts.tsv" ) -> mouse_counts

# Linke obere Ecke 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             0

Anmerkung: Falls die Tabellen in den folgenden Schritten zu groß für Ihren Computer werden, können Sie die Matrix mit folgendem Code etwas verkleinern:

mouse_counts[, c( 1:56, 100:154, 206:262 ) ] -> mouse_counts

Der folgende Code erzeugt aus den Probennamen (Spaltennamen der Count-Matrix) eine Tabelle mit allen Proben, indem er die Spaltennamen bei den Punkten aufteilt und auf drei neue Spalten verteilt. (Da manche Zeitpunkte auf “.5” enden, und an diesen Punkten nicht getrennt werden soll, müssen wir den Trenner als kryptische sog. “Regex” angeben. sep="\\.(?!5\\.)" bedeutet: Trenne an Punkten, denen kein “5.” folgt.)

tibble( sample = colnames(mouse_counts)[-1] ) %>%
separate( sample, c( "organ", "timepoint", "replicate" ), sep="\\.(?!5\\.)", remove=FALSE ) %>%
mutate( timepoint = fct_inorder(timepoint) ) -> sample_table

sample_table
# A tibble: 316 × 4
   sample        organ 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        
 7 Brain.e11.5.4 Brain e11.5     4        
 8 Brain.e12.5.1 Brain e12.5     1        
 9 Brain.e12.5.2 Brain e12.5     2        
10 Brain.e12.5.3 Brain e12.5     3        
# ℹ 306 more rows

Das Format für die Zeitpunkte der Probennahme ist wie folgt zu verstehen: e13.5 bedeutet “embryonic day 13.5”, d.h. 13.5 Tage nach Befruchtung der Eizelle und P3 bedeutet: 3 Tage seit Geburt. Das fct_inorder stellt sicher, dass die Zeitpunkte später beim Plotten in derselben Reihenfolge angeordnet werden wie in den Spalten der Matrix.

Schritt 1: Pivotieren

Wandeln Sie die Count-Matrix in eine lange Tabelle um, mit einer Zeile pro Wert und Spalten “sample”, “gene_id” und “count” für die Proben-Bezeichnung (Spaltenname in der Count-Matrix), Gen-ID und Read-Count.

Schritt 2: Normalisieren

Berechnen Sie für jede Probe die Gesamtzahl der Reads, indem Sie nach Probe (“sample”) gruppieren und dann summieren. Teilen Sie die Readzahl für jedes Gen durch die Gesamzahl und multiplizieren Sie mit 1 Million. Nennen Sie die Spalte mit diesen Werten “cpm”, für “counts per million”.

Tipp: Sie brauchen kein summerize. Verwenden Sie group_by und danach ein mutate mit count / sum(count).

Schritt 3: Table-Joins

Ergänzen Sie dann die Spalten für Organ, Zeitpunkt und Replikat, indem Sie die Probentabelle per Table-Join anfügen.

Hier finden Sie eine Tabelle, die die Gen-IDs den Gen-Symbolen (Gen-Namen) zuordnet. Fügen Sie auch diese Tabelle hinzu.

Schritt 4: Plotten

Plotten für ein Gen die Expression (CPM-Werte, mit logarithmischer Achse) gegen die Zeitpunkte, facettiert nach Organ. Unten sehen Sie zum Beispiel der Plot für das Gen “Alb”, das für das Protein Serum-Albumin kodiert.

Versuchen Sie folgende Gene:

  • Alb
  • Hbb-b1 und Hbb-bh3
  • Mbp
  • Islr2
  • Gpt

Schlagen Sie jeweils nach, was über die Gene bekannt ist.

gene <- "Alb"

long_tbl %>%
filter( gene_name == gene ) %>%
ggplot() +
  geom_point( aes( x=timepoint, y=cpm ), size=.3 ) +
  facet_wrap( ~ organ ) +
  scale_y_log10() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  ggtitle( gene )