Einführung in R Programming

`Dr. Dr. Andreas Mock`

Nationales Centrum für Tumorerkrankungen (NCT) Heidelberg

Kursinformationen

Wintersemester 2021/2022 jeweils 15:00-16:30 Uhr

Termin 1 (Intro): 12.10.21 | INF 400 SR2
Termin 2 (Transformation & Plotting): 19.10.21 | INF 400 SR2
Termin 3 (Input & Statistik): 23.11.21 | INF 400 SR1
Termin 4 (Überlebenszeitanalyse): 14.12.21 | INF 400 SR2
Termin 5 (RNA-seq Analyse): 08.02.21 | INF 400 SR1
Termin 6 (scRNA-seq Analyse): 22.02.21 | INF 400 SR2
Termin 7 (Bildanalyse & Machine Learning): 08.03.21 | INF 400 SR2

Vorbereitungen

Installation von R und RStudio

Zusätzliches Kursangebot - Data Science Seminars

Seminar 1 (Welche Programmiersprache als Mediziner erlernen? Pros und Cons von R und Python.) | 09.11.21 | 15-16:30 Uhr | INF 400 SR2
Seminar 2 (Weiterbildungsmöglichkeiten und Karrierewege) | 07.12.21 | 15-16:30 Uhr | INF 400 SR2

1 | R 101

Lehrkapitel

Präsentation (gleicher Inhalt wie Website)

Einführung

Was ist R?

R ist eine freie Programmiersprache für statistische Berechnungen und Grafiken.
Obwohl R bereits alt ist (Erscheinungsjahr 1993) gilt diese als Standardsprache für statistische Problemstellungen in Wirtschaft und Wissenschaft
16102 Formelsammlungen (Stand 11.06.2019) für spezifische Fragestellungen (sog. Pakete)
Hoch-qualitative und individuelle Grafiken - viele Wissenschaftler benutzen R nur hierzu
Sowohl R, als auch alle Pakete sind kostenlos!!

RStudio

Grafische Benutzeroberfläche und Entwicklungsumgebung für R

Abbildung: R for Data Science, Hadley Wickham & Garrett Grolemund 2016

Data Science

Data Science mit R umfasst den Import, Aufbereitung, Transformation, Visualisierung, Modellierung und Kommunikation von Daten.

Abbildung: R for Data Science, Hadley Wickham & Garrett Grolemund 2016

Objekttypen

R ist eine Objekt-orientierte Programmiersprache. Es dreht sich daher im Grunde alles darum Objekte herzustellen, zu manipulieren und zu visualisieren.

# Zuweisung von Zahlen zu einem Objekt
alter <- 67

R ist case-sensitive - Groß- und Kleinschreibung ist wichtig!

# Objekt `alter` ausgeben
alter

## [1] 67

# .. entspricht nicht
Alter

## Error in eval(expr, envir, enclos): object 'Alter' not found

Vektor

Def: Sammlung mehrerer Objekte gleicher Art (Länge 1 ist möglich).

# numerischer Vektor mit Länge 1
x1 <- 5 

# Charaktervektor mit Länge 1 (Text wird in "" gesetzt)
x2 <- "green"

# Vektoren mit Länge 3
y1 <- c(1,3,9)
y2 <- c("gene1","gene2","gene3")

Subsetting:

# um den vollständigen Vektor y2 auszugeben
y2

## [1] "gene1" "gene2" "gene3"

# um nur die ersten beiden Einträge des Vektors y2 auszugeben
y2[1:2]

## [1] "gene1" "gene2"

Matrix

Kombination mehrerer Vektoren gleichen Typs (numerisch oder Charakter). Die Matrix kann Zeilen- und Spaltennamen haben.

matrix <-  cbind(y1, y1, y1)
rownames(matrix)  <-  y2
colnames(matrix) <-  c("sample1","sample2","sample3")
matrix

##       sample1 sample2 sample3
## gene1       1       1       1
## gene2       3       3       3
## gene3       9       9       9

Ein Subset kann man sich mit der folgenden Syntax anzeigen lassen:

matrix[Zeile,Spalte]

Bespiele hierfür sind:

matrix[1,]

## sample1 sample2 sample3 
##       1       1       1

matrix[,3]

## gene1 gene2 gene3 
##     1     3     9

matrix[1:2,]

##       sample1 sample2 sample3
## gene1       1       1       1
## gene2       3       3       3

Dataframe

Im Gegensatz zu Matrizen können in Dataframes Vektoren verschiedenen Typs (z.B. numerischer Vektor und Charaktervektor) miteinander kombiniert werden. Wichtig: Die Vektoren müssen die gleiche Länge haben.

df <- data.frame(age=c(50,47,87),
                gender=c("male","male","female"))
df

##   age gender
## 1  50   male
## 2  47   male
## 3  87 female

Somit eignen sich Dataframes insbesondere für die Analyse von Patientenmetadaten im Rahmen von molekularbiologischen Experimenten oder klinischen Studien.

df

##   age gender
## 1  50   male
## 2  47   male
## 3  87 female

Eine Besonderheit von Dataframes ist die Möglichkeit einzelne Spalten durch den Spaltennamen zu selektieren.

df$age

## [1] 50 47 87

Dies entspricht der folgenden Matrixnotation

df[,1]

## [1] 50 47 87

Übersicht Objekttypen

Funktionen

Die Grundsyntax einer jeden Funktion ist

function(Objekt, Argumente)

Die Argumente sind hierbei fakultativ. R besitzt eine Vielzahl von Funktionen, ohne dass zusätzliche Packete geladen werden müssen.

sum(y1)

## [1] 13

mean(y1)

## [1] 4.333333

Die Funktion help öffnet die Dokumentation in RStudio und zeigt die notwendigen Objekte und Argumente zu jeder Funktion an. Als Beispiel, was genau macht die Funktion cbind?

help(cbind)

## starting httpd help server ... done

Objekttyp ausgeben

Die Funktion class ermöglicht es den Typ eines Objektes zu eruieren:

class(df)

## [1] "data.frame"

Pakete

Ein Paket ist nicht anderes als eine wohldurchdachte Formelsammlung, die für eine spezifische wissenschaftliche Fragestellung (z.B. die Analyse von Sequenzierungsdaten) entwickelt wurde.

Installation

install.packages("tidyverse")

Ins Environment laden

library(tidyverse)

Bespieldaten des Kurses

Metadaten des The Cancer Genome Atlas (TCGA) zur Analyse von Kopf-Hals-Tumoren (head and neck squamous cell carcinoma; HNSCC). Der Datensatz fasst die wichtigsten klinisch-pathologischen Charakteristika der Studienkohorte (n=279) zusammen.

Link zur Originalpublikation

#Datensatz in R laden
load(url("http://andreasmock.github.io/data/hnscc.RData"))

hnscc

## # A tibble: 279 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    69 YES               461 MALE   Oral Tongue   G3            51
##  2 TCGA~    39 YES               415 MALE   Larynx        G2            30
##  3 TCGA~    45 YES              1134 FEMALE Base of Tong~ G2            30
##  4 TCGA~    83 NO                276 MALE   Larynx        G2            75
##  5 TCGA~    47 YES               248 MALE   Floor of Mou~ G2            60
##  6 TCGA~    72 YES               190 MALE   Buccal Mucosa G1            20
##  7 TCGA~    56 YES               845 MALE   Alveolar Rid~ G2            NA
##  8 TCGA~    51 YES              1761 MALE   Tonsil        G2            NA
##  9 TCGA~    54 YES               186 MALE   Larynx        G2            62
## 10 TCGA~    58 YES               179 FEMALE Floor of Mou~ G3            60
## # ... with 269 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Funktionen zur Exploration von Dataframes

colnames(hnscc)

##  [1] "id"            "age"           "alcohol"       "days_to_death"
##  [5] "gender"        "neoplasm_site" "grade"         "pack_years"   
##  [9] "tabacco_group" "tumor_stage"   "vital_status"

head(hnscc$age)

## [1] 69 39 45 83 47 72

summary(hnscc$age)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.00   53.00   61.00   61.32   69.00   90.00

head(hnscc$alcohol)

## [1] "YES" "YES" "YES" "NO"  "YES" "YES"

table(hnscc$alcohol)

## 
##  NO YES 
##  85 188

table(is.na(hnscc$alcohol))

## 
## FALSE  TRUE 
##   273     6

summary(hnscc$days_to_death)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   218.8   443.0   789.0   999.2  6416.0       1

table(hnscc$gender)

## 
## FEMALE   MALE 
##     76    203

table(hnscc$neoplasm_site)

## 
## Alveolar Ridge Base of Tongue  Buccal Mucosa Floor of Mouth    Hard Palate 
##              7             12              8             26              5 
##    Hypopharynx         Larynx            Lip    Oral Cavity    Oral Tongue 
##              2             72              1             49             76 
##     Oropharynx         Tonsil 
##              2             19

table(hnscc$grade)

## 
##  G1  G2  G3  G4  GX 
##  23 176  71   1   8

summary(hnscc$pack_years)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
##   0.01685  30.00000  45.00000  50.62485  60.00000 300.00000       125

table(hnscc$tabacco_group)

## 
## Current reformed smoker for < or = 15 years 
##                                          81 
##      Current reformed smoker for > 15 years 
##                                          49 
##                              Current smoker 
##                                          90 
##                         Lifelong Non-smoker 
##                                          52

table(hnscc$tumor_stage)

## 
##   Stage I  Stage II Stage III Stage IVA Stage IVB 
##        14        44        38       139         5

table(hnscc$vital_status)

## 
## DECEASED   LIVING 
##      116      163

Zusammenfassung der Funktionen

# Tabelle von kategoriellen Daten
table(<data>)

# Tabelle von fehlenden Informationen
table(is.na(<data>))

# Summary von kontinuierlichen Daten
summary(<data>)

Literaturempfehlungen

Programmieren mit R
Hands-On Programming with R, Garrett Grolemund
- Kapitel 1, 2, 3, 4, 5, 7, 9
Link

Die “Bibel” für R User
R for Data Science, Hadley Wickham & Garrett Grolemund
Link

Kochbuch für R Plots
R Graphics Cookbock, Winston Chang
Link

Übung

Installiert das Paket ‘tidyverse’ und ladet es in euer Environment.
Ladet den Beispieldatansatz in euer Environment

load(url("http://andreasmock.github.io/data/hnscc.RData"))

Beantwortet die folgenden Fragen zum Datensatz

Wie viele Patienten trinken regelmäßig Alkohol?
Was ist die häufigste Lokalisation der Tumoren?
Bei wie vielen Patienten ist das Grading unbekannt?
Was ist die maximale Anzahl von Packyears, die angegeben wurde?
Wie viele Angaben fehlen in der Spalte “pack_years”?
In welchem Bereich bewegte sich das Überleben der Patienten?

Lösungen

Link zu Lösungen der Übung

2 | Transformation und Plotting

Lehrkapitel

Präsentation (gleicher Inhalt wie Website)

Objekttypen 2.0

Faktor

Spezielle Unterform eines Vektors für kategoriale Variablen. Ein Faktor fasst die Kategorien (= Levels) zusammen.

# Vektor
head(hnscc$grade)

## [1] "G3" "G2" "G2" "G2" "G2" "G1"

# Faktor
head(as.factor(hnscc$grade))

## [1] G3 G2 G2 G2 G2 G1
## Levels: G1 G2 G3 G4 GX

Tibbles

Tibbles sind morderne Dataframes.

class(hnscc)

## [1] "tbl_df"     "tbl"        "data.frame"

Der hnscc Datensatz ist eigentlich streng genommen kein Dataframe, sondern ein sogenannter Tibble, die moderne Weiterentwicklung eines R Dataframes.

hnscc

## # A tibble: 279 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    69 YES               461 MALE   Oral Tongue   G3            51
##  2 TCGA~    39 YES               415 MALE   Larynx        G2            30
##  3 TCGA~    45 YES              1134 FEMALE Base of Tong~ G2            30
##  4 TCGA~    83 NO                276 MALE   Larynx        G2            75
##  5 TCGA~    47 YES               248 MALE   Floor of Mou~ G2            60
##  6 TCGA~    72 YES               190 MALE   Buccal Mucosa G1            20
##  7 TCGA~    56 YES               845 MALE   Alveolar Rid~ G2            NA
##  8 TCGA~    51 YES              1761 MALE   Tonsil        G2            NA
##  9 TCGA~    54 YES               186 MALE   Larynx        G2            62
## 10 TCGA~    58 YES               179 FEMALE Floor of Mou~ G3            60
## # ... with 269 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Datentransformation

HNSCC Datensatz laden.

load(url("http://andreasmock.github.io/data/hnscc.RData"))

Warum Datentransformation? Oftmals sind wir nur an Teilmengen eines Datensatzes interessiert, bzw. möchten Proben nach verschiedenen Merkmalen zusammengruppieren.

Die Funktionen zur Datentransformation sind innerhalb des “Tidyverse” im dplyr Paket zu finden. Habt ihr das tidyverse Paket geladen, so wird automatisch auch das dplyr Paket geladen.

library(tidyverse)

Übersicht

Transformation	Funktion
Zeilen filtern	`filter()`
Zeilen sortieren	`arrange()`
Spalten selektieren	`select()`
Spaltennamen umbenennen	`rename()`
Neue Spalten hinzufügen	`mutate()`
Gruppenweise transformieren	`group_by()` & `summarize()`
Transformationen kombinieren	pipe Funktion `%>%`

Zeilen filtern mit `filter()`

Die Funktion filter() ermöglicht es uns ein Subset aus den Zeilen auszuwählen. Das erste Argument ist das Objekt, die weiteren Argumente sind die Spalten, wonach wir filtern möchten.

young <- filter(hnscc, age<50)

young

## # A tibble: 42 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    39 YES               415 MALE   Larynx        G2            30
##  2 TCGA~    45 YES              1134 FEMALE Base of Tong~ G2            30
##  3 TCGA~    47 YES               248 MALE   Floor of Mou~ G2            60
##  4 TCGA~    41 YES               242 FEMALE Oral Tongue   G2            NA
##  5 TCGA~    47 YES               395 MALE   Floor of Mou~ G2            40
##  6 TCGA~    28 YES               113 MALE   Oral Tongue   G2             1
##  7 TCGA~    48 NO               2891 MALE   Tonsil        G3            NA
##  8 TCGA~    19 NO                240 MALE   Oral Tongue   G2            NA
##  9 TCGA~    48 YES               397 FEMALE Oral Tongue   G3            20
## 10 TCGA~    48 YES               252 MALE   Larynx        G3            15
## # ... with 32 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

larynx <- filter(hnscc, neoplasm_site=="Larynx")

larynx

## # A tibble: 72 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    39 YES               415 MALE   Larynx        G2            30
##  2 TCGA~    83 NO                276 MALE   Larynx        G2            75
##  3 TCGA~    54 YES               186 MALE   Larynx        G2            62
##  4 TCGA~    53 YES               152 MALE   Larynx        G2            60
##  5 TCGA~    62 YES               244 MALE   Larynx        G2            46
##  6 TCGA~    60 YES               450 FEMALE Larynx        G2            40
##  7 TCGA~    68 YES               186 MALE   Larynx        G3            60
##  8 TCGA~    67 YES               412 MALE   Larynx        G2            NA
##  9 TCGA~    56 YES               194 MALE   Larynx        G2            80
## 10 TCGA~    52 YES               369 MALE   Larynx        G3           120
## # ... with 62 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Logische Operatoren

Die doppelten Gleichheitszeichen entsprechen der Frage: Ist der Eintrag in neoplasm_site = "Larynx". Das Resultat der Frage ist ein Vektor mit den Informationen TRUE oder FALSE pro Eintrag eines Vektors.

table(hnscc$neoplasm_site=="Larynx")

## 
## FALSE  TRUE 
##   207    72

Die Notation um alle Sites außer Larynx zu filtern ist, ein Ausrufezeichen vor den Ausdruck zu setzen:

filter(hnscc, !neoplasm_site=="Larynx")

Mehrere Sites können wie folgt ausgewählt werden:

filter(hnscc, neoplasm_site %in% c("Tonsil", "Oral Tongue", "Hard Palate"))

Im Filterprozess können Informationen aus beliebig vielen Spalten miteinander kombiniert werden.

young_larynx <- filter(hnscc, age<50, neoplasm_site=="Larynx")

young_larynx

## # A tibble: 8 x 11
##   id      age alcohol days_to_death gender neoplasm_site grade pack_years
##   <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
## 1 TCGA~    39 YES               415 MALE   Larynx        G2            30
## 2 TCGA~    48 YES               252 MALE   Larynx        G3            15
## 3 TCGA~    49 <NA>              201 MALE   Larynx        G3            NA
## 4 TCGA~    47 YES                42 MALE   Larynx        G3            40
## 5 TCGA~    45 NO                 93 FEMALE Larynx        G2            60
## 6 TCGA~    49 NO                600 MALE   Larynx        G3            16
## 7 TCGA~    38 NO                669 MALE   Larynx        GX            21
## 8 TCGA~    47 YES                35 MALE   Larynx        G2            20
## # ... with 3 more variables: tabacco_group <chr>, tumor_stage <chr>,
## #   vital_status <chr>

Zeilen sortieren mit `arrange()`

Die Funktion arrange() sortiert Zeilen nach Spalteninformationen.

arrange(hnscc, grade)

## # A tibble: 279 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    72 YES               190 MALE   Buccal Mucosa G1            20
##  2 TCGA~    65 YES              1635 MALE   Hard Palate   G1            NA
##  3 TCGA~    61 YES               236 MALE   Oral Tongue   G1            46
##  4 TCGA~    55 YES               413 FEMALE Floor of Mou~ G1            60
##  5 TCGA~    52 YES              1440 MALE   Oral Cavity   G1            45
##  6 TCGA~    45 YES               759 MALE   Oral Tongue   G1            NA
##  7 TCGA~    69 YES              1430 MALE   Oral Cavity   G1            54
##  8 TCGA~    36 YES               913 FEMALE Oral Tongue   G1            NA
##  9 TCGA~    67 NO                946 FEMALE Oral Tongue   G1            30
## 10 TCGA~    62 NO                743 MALE   Oral Tongue   G1            NA
## # ... with 269 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Hierbei kann wie auch beim Filtern eine Sortierung in mehreren Schritten erfolgen.

arrange(hnscc, age, grade)

## # A tibble: 279 x 11
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    19 NO                240 MALE   Oral Tongue   G2            NA
##  2 TCGA~    26 YES               908 MALE   Oral Tongue   G2            NA
##  3 TCGA~    26 YES              1315 MALE   Oral Tongue   G2            NA
##  4 TCGA~    28 YES               113 MALE   Oral Tongue   G2             1
##  5 TCGA~    29 <NA>              761 FEMALE Oral Tongue   GX            NA
##  6 TCGA~    32 YES                64 FEMALE Oral Tongue   G2            NA
##  7 TCGA~    34 YES               327 MALE   Oral Tongue   G2            NA
##  8 TCGA~    35 YES              1152 FEMALE Tonsil        GX            NA
##  9 TCGA~    36 YES               913 FEMALE Oral Tongue   G1            NA
## 10 TCGA~    38 YES               351 MALE   Tonsil        G2            26
## # ... with 269 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Spalten selektieren mit `select()`

select(hnscc, days_to_death, vital_status)

## # A tibble: 279 x 2
##    days_to_death vital_status
##            <int> <chr>       
##  1           461 DECEASED    
##  2           415 DECEASED    
##  3          1134 DECEASED    
##  4           276 DECEASED    
##  5           248 LIVING      
##  6           190 LIVING      
##  7           845 LIVING      
##  8          1761 DECEASED    
##  9           186 LIVING      
## 10           179 LIVING      
## # ... with 269 more rows

Umgekehrt können auch Spalten ausgeschlossen werden

select(hnscc, -c(id,age))

## # A tibble: 279 x 9
##    alcohol days_to_death gender neoplasm_site grade pack_years tabacco_group
##    <chr>           <int> <chr>  <chr>         <chr>      <dbl> <chr>        
##  1 YES               461 MALE   Oral Tongue   G3            51 Current smok~
##  2 YES               415 MALE   Larynx        G2            30 Current smok~
##  3 YES              1134 FEMALE Base of Tong~ G2            30 Current refo~
##  4 NO                276 MALE   Larynx        G2            75 Current refo~
##  5 YES               248 MALE   Floor of Mou~ G2            60 Current smok~
##  6 YES               190 MALE   Buccal Mucosa G1            20 Current refo~
##  7 YES               845 MALE   Alveolar Rid~ G2            NA Lifelong Non~
##  8 YES              1761 MALE   Tonsil        G2            NA Lifelong Non~
##  9 YES               186 MALE   Larynx        G2            62 Current refo~
## 10 YES               179 FEMALE Floor of Mou~ G3            60 Current refo~
## # ... with 269 more rows, and 2 more variables: tumor_stage <chr>,
## #   vital_status <chr>

Spaltennamen umbenennen mit `rename()`

rename(hnscc, barcode=id)

## # A tibble: 279 x 11
##    barcode   age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr>   <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA-B~    69 YES               461 MALE   Oral Tongue   G3            51
##  2 TCGA-B~    39 YES               415 MALE   Larynx        G2            30
##  3 TCGA-B~    45 YES              1134 FEMALE Base of Tong~ G2            30
##  4 TCGA-B~    83 NO                276 MALE   Larynx        G2            75
##  5 TCGA-B~    47 YES               248 MALE   Floor of Mou~ G2            60
##  6 TCGA-B~    72 YES               190 MALE   Buccal Mucosa G1            20
##  7 TCGA-B~    56 YES               845 MALE   Alveolar Rid~ G2            NA
##  8 TCGA-B~    51 YES              1761 MALE   Tonsil        G2            NA
##  9 TCGA-B~    54 YES               186 MALE   Larynx        G2            62
## 10 TCGA-B~    58 YES               179 FEMALE Floor of Mou~ G3            60
## # ... with 269 more rows, and 3 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>

Neue Spalten hinzufügen mit `mutate()`

hnscc <- mutate(hnscc, years_to_death=(days_to_death/365))

summary(hnscc$years_to_death)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.0000  0.5993  1.2137  2.1615  2.7377 17.5781       1

Gruppenweise transformieren `group_by()` und `summarize()`

by_site <- group_by(hnscc, neoplasm_site)

summarize(by_site, mean_age=mean(age))

## # A tibble: 12 x 2
##    neoplasm_site  mean_age
##    <chr>             <dbl>
##  1 Alveolar Ridge     67.7
##  2 Base of Tongue     61.7
##  3 Buccal Mucosa      70.6
##  4 Floor of Mouth     63.4
##  5 Hard Palate        76.4
##  6 Hypopharynx        59.5
##  7 Larynx             61.1
##  8 Lip                69  
##  9 Oral Cavity        66.5
## 10 Oral Tongue        56.9
## 11 Oropharynx         55.5
## 12 Tonsil             53.5

Transformationen kombinieren mit der pipe Funktion `%>%`

hnscc  %>%
    filter(!neoplasm_site %in% c("Base of Tongue","Oral Tongue")) %>%
    group_by(neoplasm_site) %>%
    summarize(count=n(),
              mean_age=mean(age))

## # A tibble: 10 x 3
##    neoplasm_site  count mean_age
##    <chr>          <int>    <dbl>
##  1 Alveolar Ridge     7     67.7
##  2 Buccal Mucosa      8     70.6
##  3 Floor of Mouth    26     63.4
##  4 Hard Palate        5     76.4
##  5 Hypopharynx        2     59.5
##  6 Larynx            72     61.1
##  7 Lip                1     69  
##  8 Oral Cavity       49     66.5
##  9 Oropharynx         2     55.5
## 10 Tonsil            19     53.5

Datenvisualisierung mit ggplot2

Funktionsweise der ggplot Funktion

Leere Leinwand. age auf der x-Achse und days_to_death auf der y-Achse.

ggplot(hnscc, aes(x=age, y=days_to_death))

ggplot(hnscc, aes(x=age, y=days_to_death))

Mit den sogenannten Aesthetics aes definieren wir die Dimensionen an Informationen, die wir im Plot darstellen möchten.

Dieser leeren Leinwand werden nun sogenannte geoms hinzugefügt, z.B. geom_point für einen Dotplot.

ggplot(hnscc, aes(x=age, y=days_to_death)) +
    geom_point()

Dotplot

ggplot(hnscc, aes(x=age, y=days_to_death)) +
    geom_point()

## Warning: Removed 1 rows containing missing values (geom_point).

Boxplot

ggplot(hnscc, aes(x=grade, y=pack_years)) +
    geom_boxplot()

## Warning: Removed 125 rows containing non-finite values (stat_boxplot).

Histogramm

ggplot(hnscc, aes(x=age)) +
    geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Density plot

ggplot(hnscc, aes(x=age)) +
    geom_density()

Transformation für Barplot

hnscc  %>%
    group_by(neoplasm_site) %>%
    summarize(mean_age=mean(age)) %>%
    ggplot(aes(x=reorder(neoplasm_site,mean_age),y=mean_age)) +
    geom_bar(stat="identity") +
    coord_flip() +
    xlab("anatomical site") +
    ylab("mean age")

Aesthetics

Bisher haben wir als aesthetics nur die x- und y-Achse verwendet. ggplot2 bietet jedoch noch weitere Dimensionen von Daten als aesthetics zu definieren

Aesthetic - color

Coloring - kategoriale Variable neoplasm site.

ggplot(hnscc, aes(x=age, y=days_to_death, color=neoplasm_site)) +
    geom_point() +
    guides(color=guide_legend(ncol=2))

## Warning: Removed 1 rows containing missing values (geom_point).

Coloring - numerische Variable packyears.

ggplot(hnscc, aes(x=age, y=days_to_death, color=pack_years)) +
    geom_point()

## Warning: Removed 1 rows containing missing values (geom_point).

Aesthetic - size

ggplot(hnscc, aes(x=age, y=days_to_death, size=pack_years)) +
    geom_point()

## Warning: Removed 126 rows containing missing values (geom_point).

Aesthetic - shape

ggplot(hnscc, aes(x=age, y=days_to_death, shape=grade)) +
    geom_point()

## Warning: Removed 1 rows containing missing values (geom_point).

Facetting

Über die Aesthetics hinaus gibt es die Möglichkeit Plots nach kategorialen Variablen zu stratefizieren.

ggplot(hnscc, aes(x=age, y=days_to_death, color=neoplasm_site)) +
    geom_point() +
    facet_wrap(~alcohol) +
    guides(color=guide_legend(ncol=2))

## Warning: Removed 1 rows containing missing values (geom_point).

Weiterführende Informationen

Datentransformation mit dplyr
Link: http://r4ds.had.co.nz/transform.html

Datenvisualisierung mit ggplot2
Link: http://r4ds.had.co.nz/data-visualisation.html

Übung

Vergleicht den Output eines “normalen” Dataframes mit einem Tibble. Welche zusätzlichen Informationen bzw. Funktionen beeinhaltet ein Tibble?

#Dataframe
head(as.data.frame(hnscc))

##             id age alcohol days_to_death gender  neoplasm_site grade pack_years
## 1 TCGA-BA-4074  69     YES           461   MALE    Oral Tongue    G3         51
## 2 TCGA-BA-4076  39     YES           415   MALE         Larynx    G2         30
## 3 TCGA-BA-4077  45     YES          1134 FEMALE Base of Tongue    G2         30
## 4 TCGA-BA-4078  83      NO           276   MALE         Larynx    G2         75
## 5 TCGA-BA-5149  47     YES           248   MALE Floor of Mouth    G2         60
## 6 TCGA-BA-5151  72     YES           190   MALE  Buccal Mucosa    G1         20
##                                 tabacco_group tumor_stage vital_status
## 1                              Current smoker   Stage IVA     DECEASED
## 2                              Current smoker        <NA>     DECEASED
## 3 Current reformed smoker for < or = 15 years   Stage IVA     DECEASED
## 4 Current reformed smoker for < or = 15 years        <NA>     DECEASED
## 5                              Current smoker   Stage IVA       LIVING
## 6      Current reformed smoker for > 15 years   Stage IVA       LIVING
##   years_to_death
## 1      1.2630137
## 2      1.1369863
## 3      3.1068493
## 4      0.7561644
## 5      0.6794521
## 6      0.5205479

#Tibble
hnscc

## # A tibble: 279 x 12
##    id      age alcohol days_to_death gender neoplasm_site grade pack_years
##    <chr> <int> <chr>           <int> <chr>  <chr>         <chr>      <dbl>
##  1 TCGA~    69 YES               461 MALE   Oral Tongue   G3            51
##  2 TCGA~    39 YES               415 MALE   Larynx        G2            30
##  3 TCGA~    45 YES              1134 FEMALE Base of Tong~ G2            30
##  4 TCGA~    83 NO                276 MALE   Larynx        G2            75
##  5 TCGA~    47 YES               248 MALE   Floor of Mou~ G2            60
##  6 TCGA~    72 YES               190 MALE   Buccal Mucosa G1            20
##  7 TCGA~    56 YES               845 MALE   Alveolar Rid~ G2            NA
##  8 TCGA~    51 YES              1761 MALE   Tonsil        G2            NA
##  9 TCGA~    54 YES               186 MALE   Larynx        G2            62
## 10 TCGA~    58 YES               179 FEMALE Floor of Mou~ G3            60
## # ... with 269 more rows, and 4 more variables: tabacco_group <chr>,
## #   tumor_stage <chr>, vital_status <chr>, years_to_death <dbl>

Erstellt die folgenden Teilmengen des Datensatzes, bzw. führt die folgenden Transformationen aus:

Junge Patienten (< 50 Jahre), Nichtraucher, Zeilen sortiert nach Alter.
Dataframe mit Spalten “age” und “grade”. Ändert die Spaltennamen ins Deutsche um.
Fügt dem hnscc Dataframe eine neue Spalte mit Namen “is_young” hinzu, die mit TRUE und FALSE kodiert, ob jemand < 50 Jahre ist.
Berechnet das mediane Gesamtüberleben nach Grading. Ist das Resultat biologisch plausibel?

In der heutigen Präsentation haben wir neben x-Achse und y-Achse weitere Aesthetics (shape, color, size), sowie das Facetting kennen gelernt. Überlegt euch eine sinnvolle Visualisierung (am besten Dotplot), die insgesamt 5 Dimensionen der Daten abbilden.a) In der heutigen Präsentation haben wir neben x-Achse und y-Achse weitere Aesthetics (shape, color, size), sowie das Facetting kennen gelernt. Überlegt euch eine sinnvolle Visualisierung (am besten Dotplot), die insgesamt 5 Dimensionen der Daten abbilden.
Versucht die Transformationen und Settings des folgenden Barplots nachzuvollziehen. Wie wird hier die Transformation mit dem Plotten verzahnt? Was macht z.B. coord_flip?

hnscc  %>%
    group_by(neoplasm_site) %>%
    summarize(mean_age=mean(age)) %>%
    ggplot(aes(x=reorder(neoplasm_site,mean_age),y=mean_age)) +
    geom_bar(stat="identity") +
    coord_flip() +
    xlab("anatomical site") +
    ylab("mean age")

“Googled” danach, wie man dem folgenden Plot eine Regressionsgrade hinzufügen kann:

ggplot(hnscc, aes(x=age, y=days_to_death)) +
    geom_point()

## Warning: Removed 1 rows containing missing values (geom_point).

Lösungen

Link zu Lösungen der Übung

3 | Input und Modelling

Lehrkapitel

Präsentation (gleicher Inhalt wie Website)

Dateipfade

Wie finde ich den absoluten Dateipfad heraus, in dem sich meine R Umgebung befindet, das so genannte working directory?

getwd()

## [1] "C:/Users/admmock/PowerFolders/git/content/AMhome"

Welche Dateien befinden sich in dem aktuellen Dateipfad?

list.files()

##  [1] "_site"           "_site.yml"       "AMhome.Rproj"    "data"           
##  [5] "index.Rmd"       "pics"            "pres_day1.Rmd"   "pres_day1_cache"
##  [9] "pres_day1_files" "pres_day2.Rmd"   "pres_day2_files" "pres_day3.Rmd"  
## [13] "pres_day3_files" "refs.Rmd"        "solution1.Rmd"   "solution2.Rmd"  
## [17] "solution3.Rmd"   "teaching.Rmd"    "teaching_cache"  "teaching_files"

Wie kann ich ausgehend vom working directory einen relativen Dateipfad angeben?

Mit dem Präfix ../ lässt sich eine Ordnereben, mit ../../ entsprechend 2 Ebenen nach oben gehen.

Der relative Pfad

list.files("../")

## [1] "AMhome"                         "Github_page"                   
## [3] "Github_paper_content"           "Zewener_Kammerkonzerte_Webpage"

entspricht damit dem absoluten Dateipfad

list.files("C:/Users/admmock/PowerFolders/git/content/")

Wie kann ich Daten aus dem Internet einlesen?

Einfach den Webpfad benutzen:

library(readxl)
path <- paste0("https://tcga-data.nci.nih.gov/docs/publications/sarc_2017/",
               "SARC_264_Fusion_Gene_Profiles.txt")

read_tsv(file = path)

Datenimport

Benötigte Packete: readr Paket im tidyverse, sowie das bisher noch nicht verwendete Paket readxl.

Funktionen zum Import nach Dateityp:

read_tsv: tab-separated file
read_csv: comma-separated file
read_xlsx bzw. read_xls: Excel spread sheet
read_delim: file mit beliebig anzugebendem delimiter, also Trennzeichen

Datenmodellierung

HNSCC Datensatz laden

load(url("http://andreasmock.github.io/data/hnscc.RData"))

T-Test

Der T-Test ist ein parametrischer (geht von normalverteilten Daten aus) Hypothesentest zum Vergleich von kontinuierlichen Daten zweier Gruppen.

ggplot(hnscc, aes(x=gender, y=age)) +
    geom_boxplot()

t.test(age ~ gender, data=hnscc)

## 
##  Welch Two Sample t-test
## 
## data:  age by gender
## t = 2.5518, df = 116.84, p-value = 0.01201
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.008323 7.998547
## sample estimates:
## mean in group FEMALE   mean in group MALE 
##             64.59211             60.08867

ANOVA (analysis of variance)

ANOVA bietet die Möglichkeit mehr als 2 Gruppen miteinander zu vergleichen.

ggplot(hnscc, aes(x=tabacco_group, y=age)) +
    geom_boxplot() +
    coord_flip()

summary(aov(age ~ tabacco_group, data=hnscc))

##                Df Sum Sq Mean Sq F value   Pr(>F)    
## tabacco_group   3   4385  1461.7   10.79 1.02e-06 ***
## Residuals     268  36294   135.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness

Korrelationstest

ggplot(hnscc, aes(x=age, y=pack_years)) +
    geom_point() +
    geom_smooth(method="lm")

## `geom_smooth()` using formula 'y ~ x'

cor.test(~ age + pack_years, data=hnscc)

## 
##  Pearson's product-moment correlation
## 
## data:  age and pack_years
## t = 1.7489, df = 152, p-value = 0.08233
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01811495  0.29211981
## sample estimates:
##       cor 
## 0.1404481

Häufige Missverständnisse über den p-Wert

Der Erfinder des p-Wertes

Sir Ronald Fisher (1890-1962)
Gonville & Caius College, Cambridge

”Personally, the writer prefers to set a low standard of significance at the 5 percent point […] A scientific fact should be regarded as experimentally established only if a properly designed experiment rarely fails to give this level of significance.“

Definition des p-Wertes

Wahrscheinlichkeit das gleiche Stichprobenergebnis oder ein noch extremeres zu erhalten, wenn die Nullhypothese wahr ist.

Algebraische Definition: P(X \(\geq\) x \(\mid\)~ \(H_0\)) wobei \(X\) eine Zufallsvariable und \(x\) der beobachte Wert in den Daten ist

Goodman, 2008

#1 | Wenn p<0.05, ist die Nullhypothese nur in 5% wahr

Dies ist die wohl des p-Wertes.

Der p-Wert wird unter der Annahme berechnet, dass die Nullhypothese zutrifft (P(Daten \(\mid\)~\(H_0\))), er kann daher nicht gleichzeitig die Wahrscheinlichkeit sein, dass die Nullhypothese zutrifft (P(\(H_0\) \(\mid\) Daten)).

Beispiel: Die Wahrscheinlichkeit drei Mal hintereinander Kopf beim Münzwurf zu erhalten ist p=0.125. Dies bedeutet jedoch nicht, dass die Wahrscheinlichkeit, dass die Münze fair ist nur 12.5% beträgt.

#2 | p>0.05 bedeutet, dass es keinen Unterschied zwischen den Gruppen gibt

Eine nicht signifikante Differenz bedeutet bloß, dass die beobachten Daten konsistent mit der Nullhypothese sind und nicht, dass die Nullhypothese wahrscheinlicher ist.

#3 | p=0.06 ist substantiell schlechter als p=0.04

Fisher hat den p-Wert als kontinuierliche Variable eingeführt um abzuschätzen, ob ein Ergebnis es Wert ist weiter untersucht zu werden. Die magische p-Wert Grenze von 0.05 ist völlig arbiträr. p-Werte von 0.04 und 0.06 sind sehr ähnliche Wahrscheinlichkeiten!

Goodman, 2008

#4 | Studien mit gleichem p-Wert zeigen eine ähnlich starke Effektgröße

Der folgende Plot zeigt, dass dies nicht zutrifft. Der gleiche p-Wert kann einen völlig anderen Effekt indizieren (Fig. B). Umgekehrt, kann es einen identischen Effekt bei unterschiedlichem p-Wert geben (Fig. A):

Goodman, 2008

#5 | p=0.05 bedeutet, dass man bei Wiederholung des Experiments in 5% ein nicht signifikantes Ergebnis erhält

Nur bei einer großen Effektgröße bzw. Power (i.e. Gruppengrößee) sind p-Werte bei Wiederholung des Experiments mit einer anderen Stichprobe reproduzierbar!

Literatur zum Thema

A Dirty Dozen: Twelve P-Value Misconceptions
Goodman, S
Semin Hematol. 2008 Jul;45(3):135-40.

The fickle P value generates irreproducible results
Halsey LG, Curran-Everett D, Vowler SL & Drummond GB
Nat Methods. 2015 Mar;12(3):179-85.

Übungen

Teil 1) Datenimport

Die anfänglichen Probleme mit dem Einlesen von Dateien lassen sich am besten durch die eigene Anwendung erkunden. Sucht euch txt, csv und xls/xlsx Files auf eurem Computer heraus und versucht diese in die R Umgebung zu laden. Denkt bitte auch an die Möglichkeit die Hilfe zu den einzelnen Funktionen aufzurufen. Am Ende der Übung sammeln wir die wichtigsten Probleme und Lösungsansätze.

Hier noch einmal die Übersicht über die Funktionen:

read_tsv: tab-separated file
read_csv: comma-separated file
read_xlsx bzw. read_xls: Excel spread sheet
read_delim: file mit beliebig anzugebendem delimiter, also Trennzeichen

Um Excel files zu laden, benötigt ihr noch das package readxl

install.packages("readxl")

library(tidyverse)
library(readxl)

Teil 2) Datenmodellierung

Datensatz: HNSCC

load(url("http://andreasmock.github.io/data/hnscc.RData"))

Benutzt einen Hypothesentest, um herauszufinden, ob es eine signifikante Assoziation zwischen dem Alter und dem Tumorstadium gibt.
Führt eine Korrelationsanalyse für Alter und days to death durch. Ist dies die beste Art und Weise, diese Assoziation zu eruieren, oder fällt euch eine alternative Möglichkeit ein?
Möchte man Verteilungen vergleichen, die nur sehr wenige Datenpunkte haben, ist es von Vorteil Rang-basierte Hypothesentest zu verwenden. Findet heraus, was die Rang-basierte Alternative zum T-Test ist und sucht nach der entsprechenden R-Funktion.

Lösungen

Link zu Lösungen der Übung

4 | Überlebenszeitanalysen

Lehrkapitel

Präsentation (gleicher Inhalt wie Website)

Kaum eine klinische Studie kommt ohne Überlebenszeitanalysen aus. Deshalb sollte die sachkundige Durchführung und Interpretation dieser zum Rüstzeug jedes Mediziners gehören.

Terminologie

Hazard Ratio (HR)
Censor
univariat
multivariat
Hazard Ratio (HR): Relative Wahrscheinlichkeit zwischen den Gruppen innerhalb eines Zeitintervalls (z.B. innerhalb eines Monats oder Jahres) ein Event zu haben (z.B. Tod oder Progress)
Censor: 1 = Event eingetroffen (z.B. Tod oder Progress), 0 = Event noch nicht eingetroffen
univariat: Es wird nur der Einfluss eines Faktors (z.B. Geschlecht) auf den Endpunkt (z.B. Gesamtüberleben) untersucht.
multivariat: Es wird nur der Einfluss mehrerer Faktoren auf den Endpunkt untersucht.

Pakete

survival: Paket für statistischen Berechnungen
survminer: Paket für Visualisierung der Kaplan-Meier-Kurven auf Basis von ggplot2
broom: Paket für Ergebnisstabellen

Datensatz

Wir werden wieder mit dem hnscc-Datensatz arbeiten.

Benötigte Daten

Für eine vergleichende Überlebenszeitanalysen benötigen wir drei Informationen:

die Zeit bis zu einem Event (z.B. Tod oder Progress)
der Censor (binäre Info ob Event eingetroffen ist, z.B. 1 = tod, 0 = lebend)
den zu untersuchenden Einflussfaktor (diskrete oder kontinuierliche Daten)

# 1. Zeit bis zum Tod
summary(hnscc$days_to_death)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0   218.8   443.0   789.0   999.2  6416.0       1

# 2. Censor
table(hnscc$vital_status)

## 
## DECEASED   LIVING 
##      116      163

# 3. den zu untersuchenden Einflussfaktor
colnames(hnscc)[-c(1,4,11)]

## [1] "age"            "alcohol"        "gender"         "neoplasm_site" 
## [5] "grade"          "pack_years"     "tabacco_group"  "tumor_stage"   
## [9] "years_to_death"

Datentransformation für Überlebenszeitanalyse

Der Censor muss ein numerischer Vektor aus 0 (=Event nicht eingetreten) und 1 (=Event eingetreten) sein.
Die Zeit bis zu einem Event (z.B. Tod) muss numerisch sein.

Diese Datentransformation entspricht der Übung 1 des heutigen Tages. Hierbei werden wir für die weitere Analyse ein eigenes Objekt hnscc_survival erstellen.

Erstellung eines survival Objektes

Für Überlebenszeitanalysen in R gilt es von so genannten survival Objekten gebrauch zu machen. Diese sind Teil des survival Pakets.

Die Syntax ist straightforward:

head(Surv(hnscc_survival$OS, hnscc_survival$Censor))

## [1]  461   415  1134   276   248+  190+

Kaplan-Meier-Kurve

Der Klassiker der Visualisierung von Überlebenszeitanalysen ist die so genannte Kaplan-Meier-Kurve.

Wichtig: Überlebenszeitanalysen können auch für kontinuierliche Einflussfaktoren durchgeführt werden. Hier macht natürlich eine Visualisierung mittels Kaplan-Meier-Kurve keinen Sinn.

Für diese müssen wir noch die Gruppenzugehörigkeit angeben, und das Objekt für den Plot mit survfit herstellen:

median(hnscc_survival$age)

## [1] 61

fit <- survfit(Surv(OS, Censor)~age>61, data=hnscc_survival)
fit

## Call: survfit(formula = Surv(OS, Censor) ~ age > 61, data = hnscc_survival)
## 
##    1 observation deleted due to missingness 
##                  n events median 0.95LCL 0.95UCL
## age > 61=FALSE 139     50   1591    1037    3314
## age > 61=TRUE  139     66    822     572    2166

ggsurvplot(fit, hnscc_survival)

Univariate Überlebenszeitanalyse

library(broom)

# Beispiel 1: Alter dichotomisiert nach > Median (= 61 Jahre)
tidy(coxph(Surv(OS, Censor)~age>61, data=hnscc_survival),
     exponentiate = TRUE)

## # A tibble: 1 x 7
##   term         estimate std.error statistic p.value conf.low conf.high
##   <chr>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 age > 61TRUE     1.54     0.189      2.27  0.0229     1.06      2.23

# Beispiel 2: Alter als kontiniuerliche Variable
tidy(coxph(Surv(OS, Censor)~age, data=hnscc_survival),
     exponentiate = TRUE)

## # A tibble: 1 x 7
##   term  estimate std.error statistic p.value conf.low conf.high
##   <chr>    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 age       1.02   0.00827      2.69 0.00718     1.01      1.04

Estimate = HR

# Bespiel 3: Tabacco Group

table(hnscc_survival$tabacco_group)

## 
## Current reformed smoker for < or = 15 years 
##                                          81 
##      Current reformed smoker for > 15 years 
##                                          49 
##                              Current smoker 
##                                          90 
##                         Lifelong Non-smoker 
##                                          52

tidy(coxph(Surv(OS, Censor)~tabacco_group, data=hnscc_survival),
     exponentiate = TRUE)

## # A tibble: 3 x 7
##   term                   estimate std.error statistic p.value conf.low conf.high
##   <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 tabacco_groupCurrent ~    0.572     0.290    -1.93   0.0542    0.324     1.01 
## 2 tabacco_groupCurrent ~    1.14      0.234     0.551  0.582     0.719     1.80 
## 3 tabacco_groupLifelong~    0.482     0.307    -2.38   0.0173    0.264     0.879

Multivariate Überlebenszeitanalyse

In der multivariaten Überlebenszeitanalyse kombinieren wir die Einflussfaktoren zusammen, die univariat signifikant waren. Wir wollen herausfinden, ob diese unabhängige Einflussfaktoren sind.

# Beispielmodel: Alter und Tabacco group
tidy(coxph(Surv(OS, Censor)~age + tabacco_group, data=hnscc_survival),
     exponentiate = TRUE)

## # A tibble: 4 x 7
##   term                   estimate std.error statistic p.value conf.low conf.high
##   <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 age                       1.02    0.00931     2.58  0.00986    1.01      1.04 
## 2 tabacco_groupCurrent ~    0.520   0.292      -2.23  0.0255     0.293     0.923
## 3 tabacco_groupCurrent ~    1.22    0.235       0.838 0.402      0.768     1.93 
## 4 tabacco_groupLifelong~    0.553   0.309      -1.92  0.0552     0.301     1.01

Übungen

Übung 1 - Datentransformation für Überlebenszeitanalyse

Erstellt ein neues Objekt, welches ihr hnscc_survival nennt und nehmt in diesem folgende Änderungen vor:

Nennt die Spalte days_to_death in OS um und stellt sicher, dass diese ein numerischer Vektor ist.
Erstellt eine neue Spalte namens Censor, in der ihr die Spalte vital_status als Faktor speichert (as.factor()).
Ändert die Level der neu erstellten Spalte Censor von “LIVING” und “DECEASED” in 0 und 1 um. Googled hierzu, wie man die Level eines Faktors ändern kann.
Anschließend gilt es die Spalte Censor als “normalen” Vektor (as.vector()) zu speichern und dann in einen numerischen Vektor (as.numeric()) zu überführen.

Übung 2 - Erstellung eines “survival” Objektes

Die Erstellung des survival Objektes ist gleichzeitig die Kontrolle, ob die Datentransformation geklappt hat. Der Output sollte wie folgt ausschauen:

Surv(hnscc_survival$OS, hnscc_survival$Censor)

##   [1]  461   415  1134   276   248+  190+  845+ 1761   186+  179+  242+ 1635+
##  [13] 1747+  152+  244+  450    63+  395   113+  236+  186+ 2891+  275+  140+
##  [25]  557+  412+  310+  142   194+  325+  166+  343+  383+  522+  369+  395 
##  [37]  240+  436   348+  810+  875+  397   279+  289+  173+  377   866+  252 
##  [49]  493   351   360   352   380   259   413+  274+  207+  228+  143+  154+
##  [61]  201+  159+   45+  134+   97+  276   201+  188+  224+   42+  177+   83+
##  [73]  119+  112   100+   46+ 1465+ 1207+  654    89   841+  896+ 1028+  317 
##  [85]  128   870+  984  1000+ 1057+  984+  403   128   455   814+ 1006+  693+
##  [97] 2562+  358+ 1662+ 1152+  798+ 1777+ 1521+ 1202  1050+ 1125+  563   514+
## [109]  183   311+  345+  354+  233+  379+  350+  479+  282  1435+ 1190+ 1440+
## [121] 1244+ 1090   104+   93   759+  889+   30+  972+   NA+ 1036+  606   796+
## [133]  521   997+ 1430+  822   392+  907+  913+  946+  908+  898+  929+  216+
## [145]  156+  180+ 1077+  910+ 1471+ 4115+  522  3930+ 3314  2318   584   546 
## [157] 2732+ 2663+ 2326+ 2567+ 1969+  545   406   540+ 1736+ 1478+ 1007+  215 
## [169] 1315+ 1246+  405+  743+  291  2740    65   295   166   624   143   665 
## [181]  804   342  3835+  602   366  1288   458   914   185  1640  2002   334 
## [193]  217   255   861    76   126  1971  3836+ 3380+  572   385   243   274 
## [205]  159    56  1591   392   662  2166   327  3496+ 1724+  144  1444+ 1095+
## [217]  953+  797+  577   520  2900   151   361  1459    64  1099+  560  1748 
## [229] 1081  6416  2717   294    14   695   761   789     2  1037  3058   452 
## [241] 4760   106   495  2570   600   217  4334+  505   194   669   316+  321+
## [253]  236+  828+  403   325+  250+  238+   77+  152+  153+ 1149+ 1058+  516+
## [265]  477+  548   427   360+  299+   35+   30+    2+    5+   25+    0+    0+
## [277]  160+  122+  113+

Erstellt mit Hilfe der survfit Funktion ein Objekt mit Namen fit für die Kaplan-Meier-Kurve, in der ihr das Geschlecht analysiert.

Übung 3 - Kaplan-Meier Kurve und univariate Überlebenszeitanalyse

Erstellt die Kaplan-Meier-Kurve zum Objekt fit, welches ihr in Übung 2 erstellt habt.
Findet heraus, wie man einen Risk-Table erstellen kann
Ist das Geschlecht ein signifikanter prognostischer Faktor?

Übung 4 - Multivariate Überlebenszeitanalyse

# Beispielmodel: Alter und Tabacco group
tidy(coxph(Surv(OS, Censor)~age + tabacco_group, data=hnscc_survival),
     exponentiate = TRUE)

## # A tibble: 4 x 7
##   term                   estimate std.error statistic p.value conf.low conf.high
##   <chr>                     <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 age                       1.02    0.00931     2.58  0.00986    1.01      1.04 
## 2 tabacco_groupCurrent ~    0.520   0.292      -2.23  0.0255     0.293     0.923
## 3 tabacco_groupCurrent ~    1.22    0.235       0.838 0.402      0.768     1.93 
## 4 tabacco_groupLifelong~    0.553   0.309      -1.92  0.0552     0.301     1.01

Interprätiert den Output. Welche Faktoren sind unabhängige Prädiktoren?

Zugabe

Gibt es weitere signifikante Prädiktoren jenseits von Age und Tabacco Group?

Lösungen

Link zu Lösungen der Übung

5 | RNA-seq Analyse

Lehrkapitel

Präsentation (gleicher Inhalt wie Website)

Kursmaterial zur RNA-seq Analyse

Das Kursmaterial des heutigen Tages entspricht in großen Teilen Kapitel 8 des hervorragenden Buchs Modern Statistics for Modern Biology von Susan Holmes und Wolfgang Huber.

Link zu Kapitel 8

Einführung

An diesem Kurstag möchten wir uns mit der Auswertung von Daten aus einer RNA Sequenzierung (RNA-seq) beschäftigen.

RNA-seq Daten sind so genannte high-throughput count Daten, in denen wir in vielen parallelen Analysen zählen, wieviele reads pro Gen oder Transkript im Sample detektiert wurden.

Sequenzierungsbibliothek: alle Nukleotide eines Samples, welche den Input für unser Sequenzierungsexperiment darstellen
Fragmente: Da der aktuelle technologische Goldstandard (Illumina Sequenzierung) nur Nukleotide von wenigen hundert Basenpaaren analysieren kann, müssen diese vor der Analyse fragmentiert werden
Read: Sequenz welche von einem Fragment erhoben wird. Ca. 150 bp bei Illumina. Bei RNA-seq und doppelsträngiger cDNA bezieht sich dies auf jeweils einen Strang

Zeitliche Abfolge: 1. Sequenzierung > 2. Counting

Benötigte Pakete

pasilla: Paket, welches die Bespieldaten enthält
DESeq2: DAS Paket für die Analyse von RNA-seq Daten
pheatmap: Paket zur Erstellung von Heatmaps (sehr gute Alternative ComplexHeatmap Paket)
broom: Paket für Ergebnisstabellen

Beispieldatensatz

Unseren Bespieldatensatz beziehen wir aus dem Paket pasilla

fn <- system.file("extdata", "pasilla_gene_counts.tsv",
                  package = "pasilla", mustWork = TRUE)

counts <- as.matrix(read.csv(fn, sep = "\t", row.names = "gene_id"))

Sequenzierungsdaten sollten als Matrix gespeichert sein und nicht als Dataframe bzw. Tibble!

dim(counts)

## [1] 14599     7

counts[1:3,]

##             untreated1 untreated2 untreated3 untreated4 treated1 treated2
## FBgn0000003          0          0          0          0        0        0
## FBgn0000008         92        161         76         70      140       88
## FBgn0000014          5          1          0          0        4        0
##             treated3
## FBgn0000003        1
## FBgn0000008       70
## FBgn0000014        0

Diese Matrix wird auch count table genannt. Wir haben Countdaten zu 14,599 Genen von 7 Samples. Hierbei handelt es sich um die “raw” Counts.

Herausforderungen der Datenanalyse

großer Wertebereich (0-Millionen)
Die nicht-negativen Count Daten (Integers) folgen keine symmetrischen Verteilung und damit keiner! Normalverteilung

plot(density(counts), xlim=c(0,10000), xlab = "Counts", main = "", bty="none")

Im Count Daten zwischen Samples zu vergleichen bedarf es einer Normalisierung. Faktoren sind hier u.a. die Sequenzierungsbibliothek, Sequenzierungstiefe, Genlänge, GC content.

Beispieldatensatz - Informationen zu Samples und zum Experiment

RNAi Knockdown Experiment des Splicingfaktors pasilla in D. melanogaster. Zwei Konditionen: treated und untreated. Den Annotationfile laden wir mit:

annotationFile <- system.file("extdata",
  "pasilla_sample_annotation.csv",
  package = "pasilla", mustWork = TRUE)
pasillaSampleAnno <- readr::read_csv(annotationFile)

pasillaSampleAnno <- mutate(pasillaSampleAnno,
condition <- factor(condition, levels = c("untreated", "treated")),
type <- factor(sub("-.*", "", type), levels = c("single", "paired")))

pasillaSampleAnno

## # A tibble: 7 x 8
##   file  condition type  `number of lane~ `total number o~ `exon counts`
##   <chr> <chr>     <chr>            <dbl> <chr>                    <dbl>
## 1 trea~ treated   sing~                5 35158667              15679615
## 2 trea~ treated   pair~                2 12242535 (x2)         15620018
## 3 trea~ treated   pair~                2 12443664 (x2)         12733865
## 4 untr~ untreated sing~                2 17812866              14924838
## 5 untr~ untreated sing~                6 34284521              20764558
## 6 untr~ untreated pair~                2 10542625 (x2)         10283129
## 7 untr~ untreated pair~                2 12214974 (x2)         11653031
## # ... with 2 more variables: `condition <- factor(condition, levels =
## #   c("untreated", "treated"))` <fct>, `type <- factor(sub("-.*", "", type),
## #   levels = c("single", "paired"))` <fct>

with(pasillaSampleAnno,
       table(condition, type))

##            type
## condition   paired-end single-read
##   treated            2           1
##   untreated          2           2

DESeq2 - DAS R Paket für die Analyse von RNA-seq Daten

Das DESeq2 gehört zu den weltweit meist genutzen Pakten zur Analyse von RNA-seq Datensätzen.

Bevor wir mit der Analyse starten können, gilt es aus der Matrix counts und dem Dataframe der die Sampleinfo beschreibt pasillaSampleAnno ein sogenanntes DESeqDataSet zu erstellen:

mt <- match(colnames(counts), sub("fb$", "", pasillaSampleAnno$file))
stopifnot(!any(is.na(mt)))

pasilla <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = pasillaSampleAnno[mt, ],
  design    = ~ condition)

pasilla

## class: DESeqDataSet 
## dim: 14599 7 
## metadata(1): version
## assays(1): counts
## rownames(14599): FBgn0000003 FBgn0000008 ... FBgn0261574 FBgn0261575
## rowData names(0):
## colnames(7): untreated1 untreated2 ... treated2 treated3
## colData names(8): file condition ... condition <- factor(condition,
##   levels = c("untreated", "treated")) type <- factor(sub("-.*", "",
##   type), levels = c("single", "paired"))

Durchführung der differentiellen Expressionsanalyse

Ziel: Gene identifizieren, die zwischen treated und untreated Samples unterschiedlich sind. Wir wenden hier einen Hypothesentest an, der mit dem t-Test mathematisch verwandt ist.

pasilla <- DESeq(pasilla)

pasilla

## class: DESeqDataSet 
## dim: 14599 7 
## metadata(1): version
## assays(4): counts mu H cooks
## rownames(14599): FBgn0000003 FBgn0000008 ... FBgn0261574 FBgn0261575
## rowData names(22): baseMean baseVar ... deviance maxCooks
## colnames(7): untreated1 untreated2 ... treated2 treated3
## colData names(9): file condition ... type <- factor(sub("-.*", "",
##   type), levels = c("single", "paired")) sizeFactor

Results

res <- results(pasilla)
res[order(res$padj), ] %>% head

## log2 fold change (MLE): condition untreated vs treated 
## Wald test p-value: condition untreated vs treated 
## DataFrame with 6 rows and 6 columns
##                     baseMean    log2FoldChange              lfcSE
##                    <numeric>         <numeric>          <numeric>
## FBgn0039155 730.595806139728    4.619013535193  0.168706755496974
## FBgn0025111 1501.41051323996 -2.89986418283187  0.126920467532813
## FBgn0029167 3706.11653071978  2.19700028710382 0.0969888689039936
## FBgn0003360 4343.03539692487  3.17967242210135  0.143526226262858
## FBgn0035085 638.232608936723  2.56041205175815  0.137295199514033
## FBgn0039827 261.916235943103  4.16251610026134  0.232588756794067
##                          stat                pvalue                  padj
##                     <numeric>             <numeric>             <numeric>
## FBgn0039155  27.3789482915866 4.88534640286228e-165 4.06607381110228e-161
## FBgn0025111 -22.8478844996546 1.53390416993678e-115 6.38334220319191e-112
## FBgn0029167  22.6520869036895 1.33057630038134e-113 3.69146218269131e-110
## FBgn0003360  22.1539470861445 9.55682092175738e-109 1.98853551329467e-105
## FBgn0035085  18.6489553955341  1.28768628708899e-77  2.14348259348833e-74
## FBgn0039827   17.896463086334  1.25652932702639e-71  1.74301559814011e-68

Im folgenden möchten wir die Results der differentiellen Expressionsanalyse mithilfe von Visualisierungen explorieren.

Histogramm der p-Werte

ggplot(as(res, "data.frame"), aes(x = pvalue)) +
  geom_histogram(binwidth = 0.01, fill = "Royalblue", boundary = 0) +
    theme_classic()

ggplot(as(res, "data.frame"), aes(x = pvalue)) +
  geom_histogram(binwidth = 0.005, fill = "Royalblue", boundary = 0) +
    xlim(c(0,0.1)) +
    theme_classic()

MA plot

Fold change (M-value) vs. average counts (A-value)

plotMA(pasilla, ylim = c( -2, 2))

Roten Punkte = adjustierter P-Wert < 0.1

PCA plot

Die Principle Component Analysis (PCA) ist eine Methodik zur Dimensionsreduktion. Hierbei werden für jeden Sample die Infos zu allen 14,599 Genen auf zwei Dimensionen (x- und y-Achse) reduziert. Ein Sample ist ein Punkt auf dem Plot.

pas_rlog <- rlogTransformation(pasilla)
plotPCA(pas_rlog, intgroup=c("condition", "type")) + 
    coord_fixed() +
    theme_classic()

Heatmap

Eine Heatmap aller Gene darzustellen macht keinen Sinn. Es gilt Filterkriterien zu definieren.

Filter 1) Top 30 most expressed genes.

select <- order(rowMeans(assay(pas_rlog)), decreasing = TRUE)[1:30]

pheatmap( assay(pas_rlog)[select, ],
     scale = "row",
     annotation_col = as.data.frame(
        colData(pas_rlog)[, c("condition", "type")] ))

Filter 2) Adjusted p-Value < 0.01 und absolute log2 Fold Change > 2

select2 <- (res$padj<0.01 & abs(res$log2FoldChange)>2)*1

pheatmap( assay(pas_rlog)[select2 %in% 1, ],
     scale = "row",
     annotation_col = as.data.frame(
        colData(pas_rlog)[, c("condition", "type")] ))

Export der Results als csv File

write.csv(as.data.frame(res), file = "treated_vs_untreated.csv")

Übung

Die Übung des heutigen Tages besteht im Ausführen und Nachvollziehen des R Codes der Slides.

Einführung in R Programming

Dr. Dr. Andreas Mock

Kursinformationen

1 | R 101

Lehrkapitel

Einführung

Was ist R?

RStudio

Data Science

Objekttypen

Vektor

Matrix

Dataframe

Übersicht Objekttypen

Funktionen

Objekttyp ausgeben

Pakete

Bespieldaten des Kurses

Funktionen zur Exploration von Dataframes

Literaturempfehlungen

Übung

Lösungen

2 | Transformation und Plotting

Lehrkapitel

Objekttypen 2.0

Faktor

Tibbles

Datentransformation

Übersicht

Zeilen filtern mit filter()

Zeilen sortieren mit arrange()

Spalten selektieren mit select()

Spaltennamen umbenennen mit rename()

Neue Spalten hinzufügen mit mutate()

Gruppenweise transformieren group_by() und summarize()

Transformationen kombinieren mit der pipe Funktion %>%

Datenvisualisierung mit ggplot2

Aesthetics

Facetting

Weiterführende Informationen

Übung

Lösungen

3 | Input und Modelling

Lehrkapitel

Dateipfade

Datenimport

Datenmodellierung

T-Test

ANOVA (analysis of variance)

Korrelationstest

Häufige Missverständnisse über den p-Wert

Der Erfinder des p-Wertes

Definition des p-Wertes

#1 | Wenn p<0.05, ist die Nullhypothese nur in 5% wahr

#2 | p>0.05 bedeutet, dass es keinen Unterschied zwischen den Gruppen gibt

#3 | p=0.06 ist substantiell schlechter als p=0.04

#4 | Studien mit gleichem p-Wert zeigen eine ähnlich starke Effektgröße

#5 | p=0.05 bedeutet, dass man bei Wiederholung des Experiments in 5% ein nicht signifikantes Ergebnis erhält

Literatur zum Thema

Übungen

Teil 1) Datenimport

Teil 2) Datenmodellierung

Lösungen

4 | Überlebenszeitanalysen

Lehrkapitel

Terminologie

Pakete

Datensatz

Benötigte Daten

Datentransformation für Überlebenszeitanalyse

Erstellung eines survival Objektes

Kaplan-Meier-Kurve

Univariate Überlebenszeitanalyse

Multivariate Überlebenszeitanalyse

Übungen

Übung 1 - Datentransformation für Überlebenszeitanalyse

Übung 2 - Erstellung eines “survival” Objektes

Übung 3 - Kaplan-Meier Kurve und univariate Überlebenszeitanalyse

Übung 4 - Multivariate Überlebenszeitanalyse

Zugabe

`Dr. Dr. Andreas Mock`

Zeilen filtern mit `filter()`

Zeilen sortieren mit `arrange()`

Spalten selektieren mit `select()`

Spaltennamen umbenennen mit `rename()`

Neue Spalten hinzufügen mit `mutate()`

Gruppenweise transformieren `group_by()` und `summarize()`

Transformationen kombinieren mit der pipe Funktion `%>%`