Lab 5: Data Wrangling II
Package(s)
Schedule
- 08.00 - 08.30: Recap of Lab 4
- 08.30 - 08.35: Lecture
- 08.35 - 08.45: Break
- 08.45 - 12.00: Exercises
Learning Materials
Please prepare the following materials
- Book: R4DS2e: Chapter 5 Data tidying
- Book: R4DS2e: Chapter 14 Strings
- Book: R4DS2e: Chapter 16 Factors
- Book: Chapter 19 Joins
- Video: Tidy Data and tidyr - NB! Start at 7:45 and please note:
gather()
is nowpivot_longer()
andspread()
is nowpivot_wider()
- Video: Working with Two Datasets: Binds, Set Operations, and Joins
- Video: stringr (Playlist with 7 short videos)
Unless explicitly stated, do not do the per-chapter exercises in the R4DS2e book
Learning Objectives
A student who has met the objectives of the session will be able to:
- Understand and apply the various
str_*()
functions for string manipulation - Understand and apply the family of
*_join()
functions for combining data sets - Understand and apply
pivot_wider()
andpivot_longer()
- Use factors in context with plotting categorical data using
ggplot
Exercises
Prologue
Today will not be easy! But please try to remember Hadley’s words of advice:
- “The bad news is, whenever you’re learning a new tool, for a long time, you’re going to suck! It’s gonna be very frustrating! But the good news is that that is typical and something that happens to everyone and it’s only temporary! Unfortunately, there is no way to going from knowing nothing about the subject to knowing something about a subject and being an expert in it without going through a period of great frustration and much suckiness! Keep pushing through!” - H. Wickham (dplyr tutorial at useR 2014, 4:10 - 4:48)
Intro
We are upping the game here, so expect to get stuck at some of the questions. Remember - Discuss with your group how to solve the task, revisit the materials you prepared for today and naturally, the TAs and I are happy to nudge you in the right direction. Finally, remember… Have fun!
Remember what you have worked on so far:
- RStudio
- Quarto
ggplot
filter
arrange
select
mutate
group_by
summarise
- The pipe and creating pipelines
stringr
- joining data
- pivoting data
That’s quite a lot! Well done - You’ve come quite far already! Remember to think about the above tools in the following as we will synthesise your learnings so far into an analysis!
Background
In the early 20s, the world was hit by the coronavirus disease 2019 (COVID-19) pandemic. The pandemic was caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In Denmark, the virus first confirmed case was on 27 February 2020.
While initially very little was known about the SARS-CoV-2 virus, we did know the general pathology of vira. Briefly, the virus invades the cells and hijacks the intra-cellular machinery. Using the hijacked machinery, components for new virus particles are produced, eventually being packed into the viral envelope and released from the infected cell. Some of these components, viral proteins, is broken down into smaller fragments called peptides by the proteasome. These peptides are transported into the endoplasmic reticulum by the Transporter Associated with antigen Processing (TAP) protein complex. Here, they are aided by chaperones bound to the Major Histocompatilibty Complex class I (MHC-I) and then across the Golgi apparatus they finally get displayed on the surface of the cells. Note, in humans, MHC is also called Human Leukocyte Antigen (HLA) and represents the most diverse genes. Each of us have a total of 6 HLA-alleles, 3 from the maternal and 3 from the paternal side. These are further divided into 3 classes HLA-A, HLA-B and HLA-C and the combination of these constitute the HLA-haplotype for an individual. Once the peptide is bound to the MHC class I at the cell surface and exposed, the MHC-I peptide complex can be recognised by CD8+ Cytotoxic T-Lymphocytes (CTLs) via the T-cell Receptor (TCR). If a cell displays peptides of viral origin, the CTL gets activated and via a cascade induces apoptosis (programmed cell death) of the infected cell. The process is summarised in the figure below (McCarthy and Weinberg 2015).
The data we will be working with today contains data on sequenced T-cell receptors, viral antigens, HLA-haplotypes and clinical meta data for a cohort:
- “A large-scale database of T-cell receptor beta (TCR\(\beta\)) sequences and binding associations from natural and synthetic exposure to SARS-CoV-2” (Nolan et al. 2020).
Your Task Today
Today, we will emulate the situation, where you are working as a Bioinformatician / Bio Data Scientist and you have been given the data and the task of answering these two burning questions:
- What characterises the peptides binding to the HLAs?
- What characterises T-cell Receptors binding to the pMHC-complexes?
GROUP ASSIGNMENT: Today, your assignment will be to create a micro-report on these 2 questions! (Important, see: how to)
Getting Started
First, make sure to read and discuss the feedback you got from last week’s assignment!
- Then, once again go to the R for Bio Data Science RStudio Cloud Server
- Make sure you are in your
r_for_bio_data_science
project, you can verify this in the upper right corner - In the same place as your
r_for_bio_data_science.Rproj
file and existingdata
folder, create a new folder and name itdoc
- Go to the aforementioned manuscript. Download the PDF and upload it to your new
doc
folder - Open the PDF and find the link to the data
- Go to the data site (Note, you may have to create and account to download, shouldn’t take too long) . Find and download the file
ImmuneCODE-MIRA-Release002.1.zip
(CAREFUL, do not download the superseded files) - Unpack the downloaded file
- Find the files
peptide-detail-ci.csv
andsubject-metadata.csv
and compress to.zip
files - Upload the compressed
peptide-detail-ci.csv.zip
andsubject-metadata.csv.zip
files to yourdata
folder in your RStudio Cloud session - Finally, once again, create a new Quarto document for today’s exercises, containing the sections:
- Background
- Aim
- Load Libraries
- Load Data
- Data Description
- Analysis
Creating the Micro-Report
Background
Feel free to copy paste the one stated in the background-section above
Aim
State the aim of the micro-report, i.e. what are the questions you are addressing?
Load Libraries
Load the libraries needed
Load Data
Read the two data sets into variables peptide_data
and meta_data
.
Click here for hint
Think about which Tidyverse package deals with reading data and what are the file types we want to read here?Data Description
It is customary to include a description of the data, helping the reader if the report, i.e. your stakeholder, to get an easy overview
The Subject Meta Data
Let’s take a look at the meta data:
|>
meta_data sample_n(10)
# A tibble: 10 × 30
Experiment Subject `Cell Type` `Target Type` Cohort Age Gender Race
<chr> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr>
1 eMR22 1565927 PBMC C19_cI COVID-19-Con… 65 M <NA>
2 eXL32 4423 naive_CD8 C19_cI Healthy (No … 37 F White
3 eEE240 20795 naive_CD8 C19_cI Healthy (No … 23 M White
4 eLH48 535 PBMC C19_cI COVID-19-Con… 28 M White
5 eHO131 3282 PBMC C19_cI COVID-19-Con… 58 F <NA>
6 eNL192 1524 B-CD8-_PBMC C19_cII COVID-19-Con… NA <NA> <NA>
7 eHO130 9541 PBMC C19_cI Healthy (No … 28 F White
8 eLH43 3565 PBMC C19_cI COVID-19-Con… 57 M <NA>
9 eAV100 1995 PBMC C19_cII COVID-19-Con… 29 F <NA>
10 eQD115 2513 PBMC C19_cI COVID-19-Con… 48 M <NA>
# ℹ 22 more variables: `HLA-A...9` <chr>, `HLA-A...10` <chr>,
# `HLA-B...11` <chr>, `HLA-B...12` <chr>, `HLA-C...13` <chr>,
# `HLA-C...14` <chr>, DPA1...15 <chr>, DPA1...16 <chr>, DPB1...17 <chr>,
# DPB1...18 <chr>, DQA1...19 <chr>, DQA1...20 <chr>, DQB1...21 <chr>,
# DQB1...22 <chr>, DRB1...23 <chr>, DRB1...24 <chr>, DRB3...25 <chr>,
# DRB3...26 <chr>, DRB4...27 <chr>, DRB4...28 <chr>, DRB5...29 <chr>,
# DRB5...30 <chr>
Q1: How many observations of how many variables are in the data?
Q2: Are there groupings in the variables, i.e. do certain variables “go together” somehow?
T1: Re-create this plot
Read this first:
- Think about: What is on the x-axis? What is on the y-axis? And also, it looks like we need to do some
count
ing stratified byCohort
andGender
. Recall, that we can stick together adplyr
pipeline with a call toggplot
.
Does your plot look different somehow? Consider peeking at the hint…
Click here for hint
Perhaps not everyone agrees on how to denoteNA
s in data. I have seen -99
, -11
, _
and so on… Perhaps this can be dealt with in the instance we read the data from the file? I.e. in the actual function call to your read_csv()
function. Recall, how can we get information on the parameters of a ?function
- T2: Re-create this plot
Click here for hint
Perhaps there is a function, which cancut
continuous observations into a set of bins?
STOP! Make sure you handled how NA
s are denoted in the data before proceeding, see hint below T1
- T3: Look at the data and create yet another plot as you see fit. Also skip the redundant variables
Subject
,Cell Type
andTarget Type
|>
meta_data sample_n(10)
# A tibble: 10 × 27
Experiment Cohort Age Gender Race `HLA-A...9` `HLA-A...10` `HLA-B...11`
<chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr>
1 eOX49 Healthy … 21 M White A*02:01 A*26:01 B*44:02
2 eLH57 COVID-19… NA <NA> <NA> A*01:01:01 A*03:01:01 B*07:02:01
3 eLH43 COVID-19… 57 M <NA> A*02:13 A*23:01:01 B*44:02:01
4 eQD132 COVID-19… NA <NA> <NA> A*02:03:01 A*11:01:01 B*39:01:01
5 eQD139 COVID-19… NA <NA> <NA> A*01:01:01 A*29:02:01 B*56:01:01
6 eEE217 Healthy … 32 F White A*02:01 A*02:01 B*15:01
7 eLH48 COVID-19… 28 M White A*03:01:01 A*24:02:01 B*08:01:01
8 eLH41 COVID-19… 71 F <NA> A*02:01:01 A*03:01:01 B*13:02:01
9 eQD127 COVID-19… 61 F <NA> A*02:01:01 A*03:01:01 B*27:02:01
10 eJL157 COVID-19… 27 F <NA> A*11:01:01 A*11:01:01 B*07:02:01
# ℹ 19 more variables: `HLA-B...12` <chr>, `HLA-C...13` <chr>,
# `HLA-C...14` <chr>, DPA1...15 <chr>, DPA1...16 <chr>, DPB1...17 <chr>,
# DPB1...18 <chr>, DQA1...19 <chr>, DQA1...20 <chr>, DQB1...21 <chr>,
# DQB1...22 <chr>, DRB1...23 <chr>, DRB1...24 <chr>, DRB3...25 <chr>,
# DRB3...26 <chr>, DRB4...27 <chr>, DRB4...28 <chr>, DRB5...29 <chr>,
# DRB5...30 <chr>
Now, a classic way of describing a cohort, i.e. the group of subjects used for the study, is the so-called table1
and while we could build this ourselves, this one time, in the interest of exercise focus and time, we are going to “cheat” and use an R-package, like so:
NB!: This may look a bit odd initially, but if you render your document, you should be all good!
library("table1") # <= Yes, this should normally go at the beginning!
|>
meta_data mutate(Gender = factor(Gender),
Cohort = factor(Cohort)) |>
table1(x = formula(~ Gender + Age + Race | Cohort),
data = _)
COVID-19-Acute (N=4) |
COVID-19-B-Non-Acute (N=8) |
COVID-19-Convalescent (N=90) |
COVID-19-Exposed (N=3) |
Healthy (No known exposure) (N=39) |
Overall (N=144) |
|
---|---|---|---|---|---|---|
Gender | ||||||
F | 1 (25.0%) | 4 (50.0%) | 33 (36.7%) | 1 (33.3%) | 17 (43.6%) | 56 (38.9%) |
M | 2 (50.0%) | 3 (37.5%) | 36 (40.0%) | 0 (0%) | 21 (53.8%) | 62 (43.1%) |
Missing | 1 (25.0%) | 1 (12.5%) | 21 (23.3%) | 2 (66.7%) | 1 (2.6%) | 26 (18.1%) |
Age | ||||||
Mean (SD) | 50.7 (17.0) | 43.7 (7.74) | 51.5 (15.3) | 35.0 (NA) | 33.3 (9.93) | 44.9 (15.7) |
Median [Min, Max] | 52.0 [33.0, 67.0] | 42.0 [33.0, 53.0] | 53.0 [21.0, 79.0] | 35.0 [35.0, 35.0] | 31.0 [21.0, 62.0] | 42.0 [21.0, 79.0] |
Missing | 1 (25.0%) | 1 (12.5%) | 21 (23.3%) | 2 (66.7%) | 0 (0%) | 25 (17.4%) |
Race | ||||||
African American | 1 (25.0%) | 0 (0%) | 0 (0%) | 0 (0%) | 1 (2.6%) | 2 (1.4%) |
White | 2 (50.0%) | 7 (87.5%) | 13 (14.4%) | 0 (0%) | 28 (71.8%) | 50 (34.7%) |
Asian | 0 (0%) | 0 (0%) | 3 (3.3%) | 0 (0%) | 2 (5.1%) | 5 (3.5%) |
Hispanic or Latino/a | 0 (0%) | 0 (0%) | 1 (1.1%) | 0 (0%) | 0 (0%) | 1 (0.7%) |
Native Hawaiian or Other Pacific Islander | 0 (0%) | 0 (0%) | 0 (0%) | 1 (33.3%) | 0 (0%) | 1 (0.7%) |
Black or African American | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 3 (7.7%) | 3 (2.1%) |
Mixed Race | 0 (0%) | 0 (0%) | 0 (0%) | 0 (0%) | 1 (2.6%) | 1 (0.7%) |
Missing | 1 (25.0%) | 1 (12.5%) | 73 (81.1%) | 2 (66.7%) | 4 (10.3%) | 81 (56.3%) |
Note how good this looks! If you have ever done a “Table 1” before, you know how painful they can be and especially if something changes in your cohort - Dynamic reporting to the rescue!
Lastly, before we proceed, the meta_data
contains HLA data for both class I and class II (see background), but here we are only interested in class I, recall these are denoted HLA-A
, HLA-B
and HLA-C
, so make sure to remove any non-class I, i.e. the one after, denoted D
-something.
- T4: Create a new version of the
meta_data
, which with respect to allele-data only contains information on class I and also fix the odd naming, e.g.HLA-A...9
becomesA1
oandHLA-A...10
becomesA2
and so on forB1
,B2
,C1
andC2
(Think: How can werename
variables? And here, just do it “manually” per variable). Remember to assign this new data to the samemeta_data
variable
Click here for hint
Whichtidyverse
function subsets variables? Perhaps there is a function, which somehow matches
a set of variables? And perhaps for the initiated this is compatible with regular expressions (If you don’t know what this means - No worries! If you do, see if you utilise this to simplify your variable selection)
Before we proceed, this is the data we will carry on with:
|>
meta_data sample_n(10)
# A tibble: 10 × 11
Experiment Cohort Age Gender Race A1 A2 B1 B2 C1 C2
<chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 eMR12 COVID-19-C… NA <NA> <NA> A*01… A*02… B*40… B*52… C*03… C*16…
2 eLH59 COVID-19-C… NA <NA> <NA> A*01… A*02… B*40… B*52… C*03… C*16…
3 eDH113 Healthy (N… 56 <NA> <NA> A*02… A*29… B*18… B*44… C*07… C*16…
4 eQD119 COVID-19-C… 51 M <NA> A*01… A*01… B*08… B*37… C*06… C*07…
5 eXL36 Healthy (N… 37 F White A*01… A*02… B*15… B*40… C*03… C*04…
6 eAV105 COVID-19-C… 29 F <NA> A*02… A*68… B*07… B*40… C*03… C*07…
7 eAV100 COVID-19-C… 29 F <NA> A*02… A*68… B*07… B*40… C*03… C*07…
8 eHH169 Healthy (N… 24 F Blac… A*02… A*74… B*35… B*35… C*04… C*04…
9 ePD81 COVID-19-C… 64 M <NA> A*24… A*24… B*40… B*46… C*01… C*15…
10 eHH173 Healthy (N… 50 M White A*02… A*03… B*35… B*44… C*04… C*05…
Now, we have a beautiful tidy
dataset, recall that this entails, that each row is an observation, each column is a variable and each cell holds one value.
The Peptide Details Data
Let’s start with simply having a look see:
|>
peptide_data sample_n(10)
# A tibble: 10 × 7
`TCR BioIdentity` TCR Nucleotide Seque…¹ Experiment `ORF Coverage`
<chr> <chr> <chr> <chr>
1 CASSEAEDKTTNEKLFF+TCRBV02-0… ACAAAGCTGGAGGACTCAGCC… eQD108 ORF3a
2 CATSDSRGDPIPTGELFF+TCRBV24-… ATCCCCAACCAGACAGCTCTT… eAV93 nucleocapsid …
3 CASSWGPAGELFF+TCRBV27-01+TC… ATCCTGGAGTCGCCCAGCCCC… eOX52 ORF1ab
4 CASSLAAGYGYTF+TCRBV27-01+TC… ATCCTGGAGTCGCCCAGCCCC… eXL31 ORF1ab
5 CASSLAGVGLYF+TCRBV05-01+TCR… ATGAATGTGAGCACCTTGGAG… eXL31 nucleocapsid …
6 CASSYSQAGGYNEQFF+TCRBV06-02… TCGGCTGCTCCCTCCCAAACA… eLH53 ORF1ab
7 CSATSINPGQEYGYTF+TCRBV20-X+… AGTGCCCATCCTGAAGACAGC… eXL30 ORF10
8 CASSLDRAQETQYF+TCRBV07-09+T… ATCCAGCGCACAGAGCAGGGG… eHO132 ORF1ab
9 CASSKAGLAGGSNEQFF+TCRBV21-0… ACGGAGTCAGGGGACACAGCA… eOX49 ORF1ab
10 CASSFGTGGKDGYTF+TCRBV07-09+… CAGCGCACAGAGCAGGGGGAC… eHO141 nucleocapsid …
# ℹ abbreviated name: ¹`TCR Nucleotide Sequence`
# ℹ 3 more variables: `Amino Acids` <chr>, `Start Index in Genome` <dbl>,
# `End Index in Genome` <dbl>
- Q3: How many observations of how many variables are in the data?
This is a rather big data set, so let us start with two “tricks” to handle this, first:
- Write the data back into your
data
folder, using the filenamepeptide-detail-ci.csv.gz
, note the appending of.gz
, which is automatically recognised and results in gz-compression - Now, check in your data folder, that you have two files
peptide-detail-ci.csv
andpeptide-detail-ci.csv.gz
, delete the former - Adjust your reading-the-data-code in the “Load Data”-section, to now read in the
peptide-detail-ci.csv.gz
file
Click here for hint
Just as you canread
a file, you can of course also write
a file. Note the filetype we want to write here is csv
. If you in the console type e.g. readr::wr
and then hit the Tab key, you will see the different functions for writing different filetypes
Then:
- T5: As before, let’s immediately subset the
peptide_data
to the variables of interest:TCR BioIdentity
,Experiment
andAmino Acids
. Remember to assign this new data to the samepeptide_data
variable to avoid cluttering your environment with redundant variables. Bonus: Did you know you can click theEnvironment
pane and see which variables you have?
Once again, before we proceed, this is the data we will carry on with:
|>
peptide_data sample_n(10)
# A tibble: 10 × 3
Experiment `TCR BioIdentity` `Amino Acids`
<chr> <chr> <chr>
1 eAV88 CASSHRGGNQPQHF+TCRBV13-01+TCRBJ01-05 FLWLLWPVT,FLWLLWPVTL,LWLL…
2 eQD125 CASSTRSGLDEQYF+TCRBV27-01+TCRBJ02-07 HTTDPSFLGRY
3 eXL30 CASSLAGGLNEQFF+TCRBV06-05+TCRBJ02-01 AFLLFLVLI,FLAFLLFLV,FYLCF…
4 eXL31 CASSPPPGTGTADTQYF+TCRBV13-01+TCRBJ02-03 YIFFASFYY
5 eQD110 CASSFGTGNEQFF+TCRBV13-01+TCRBJ02-01 HTTDPSFLGRY
6 eXL30 CASSYTLTEMTEAFF+TCRBV06-05+TCRBJ01-01 AFLLFLVLI,FLAFLLFLV,FYLCF…
7 eEE224 CASQGQGDNEQFF+TCRBV28-01+TCRBJ02-01 FLWLLWPVT,FLWLLWPVTL,LWLL…
8 eEE240 CASSSSRAQPQHF+TCRBV06-06+TCRBJ01-05 RQLLFVVEV
9 eOX46 CASSQTLDPYEQYF+TCRBV06-05+TCRBJ02-07 FLWLLWPVT,FLWLLWPVTL,LWLL…
10 eLH48 CASSYSSNEQFF+TCRBV06-05+TCRBJ02-01 TSDLATNNLVVMAY
Q4: Is this tidy data? Why/why not?
T6: See if you can find a way to create the below data, from the above
|>
peptide_data sample_n(size = 10)
# A tibble: 10 × 5
Experiment CDR3b V_gene J_gene `Amino Acids`
<chr> <chr> <chr> <chr> <chr>
1 eEE243 CASSPTVVVAGELFF TCRBV09-01 TCRBJ02-02 FLNGSCGSV
2 eXL27 CASSYEGGRSEQYF TCRBV27-01 TCRBJ02-07 SELVIGAVI,SELVIGAVIL
3 eXL37 CASSDTDTQYF TCRBV09-01 TCRBJ02-03 YLNTLTLAV
4 eXL30 CASSLGEMVVTEAFF TCRBV07-09 TCRBJ01-01 ALALLLLDR,GDAALALLLL,…
5 eAV91 CSARGFTTDTQYF TCRBV20-X TCRBJ02-03 MPASWVMRI
6 eOX49 CSSGQGVYEQYF TCRBV20-X TCRBJ02-07 TMADLVYAL,YTMADLVYA
7 eHO141 CASRMGQNSNQPQHF TCRBV06-05 TCRBJ01-05 AYKTFPPTEPK,KTFPPTEPK
8 eAV88 CASSQTRDRIYEQYF TCRBV04-01 TCRBJ02-07 SVLNDILSR,VLNDILSRL
9 eLH47 CASSFAEKLFF TCRBV12-03/12-04 TCRBJ01-04 AFPFTIYSL,GYINVFAFPF,…
10 eOX54 CASSLGWGGTEAFF TCRBV05-04 TCRBJ01-01 VLWAHGFEL
Click here for hint
First: Compare the two datasets and identify what happened? Did any variables “disappear” and did any “appear”? Ok, so this is a bit tricky, but perhaps there is a function toseparate
a composite (untidy) col
umn into
a set of new variables based on a sep
arator? But what is a sep
arator? Just like when you read a file with C
omma S
eparated V
alues, a separator denotes how a composite string is divided into fields. So, look for such a repeated value, which seem to indeed separate such fields. Also, be aware, that character, which can mean more than one thing, may need to be “escaped” using an initial two backslashed, i.e. “\x”, where x denotes the character needing to be “escaped”
- T7: Add a variable, which counts how many peptides are in each observation of
Amino Acids
Click here for hint
We have been working with thestringr
package, perhaps the contains a function to somehow count the number of occurrences of a given character in a string? Again, remember you can type e.g. stringr::str_
and then hit the Tab key to see relevant functions
|>
peptide_data sample_n(size = 10)
# A tibble: 10 × 6
Experiment CDR3b V_gene J_gene `Amino Acids` n_peptides
<chr> <chr> <chr> <chr> <chr> <dbl>
1 eQD134 CSGRQGADIQYF TCRBV29-01 TCRBJ02-04 ELIRQGTDY,QELI… 3
2 ePD83 CASSPGVGLSYEQYF TCRBV19-01 TCRBJ02-07 SEHDYQIGGYTEKW… 3
3 ePD85 CASSLGQGTTYEQYF TCRBV19-01 TCRBJ02-07 SEHDYQIGGYTEKW… 3
4 eOX49 CASTTIVGQGYPYEQYF TCRBV28-01 TCRBJ02-07 FLWLLWPVT,FLWL… 7
5 eQD137 CASSHHRGWDNSPLHF TCRBV25-01 TCRBJ01-06 SYFTSDYYQL,VLH… 3
6 eEE226 CSARDASHYGYTF TCRBV20-X TCRBJ01-02 VLWAHGFEL 1
7 eXL30 CASSPGTVSNQPQHF TCRBV02-01 TCRBJ01-05 DFLEYHDVR,EDFL… 5
8 eEE240 CASSIDPTAKNIQYF TCRBV19-01 TCRBJ02-04 QYIKWPWYI,YEQY… 3
9 eEE226 CASDRGREKLFF TCRBV06-05 TCRBJ01-04 FLNGSCGSV 1
10 eHO136 CSASPGGNLYEQYF TCRBV20-X TCRBJ02-07 LLTDEMIAQY,LTD… 4
- T8: Re-create the following plot
Q4: What is the maximum number of peptides assigned to one observation?
T9: Using the
str_c()
and theseq()
functions, re-create the below
[1] "peptide_1" "peptide_2" "peptide_3" "peptide_4" "peptide_5"
Click here for hint
If you’re uncertain on how a function works, try going into the console and in this case e.g. typestr_c("a", "b")
and seq(from = 1, to = 3)
and see if you combine these?
- T10: Use, what you learned about separating in T6 and the vector-of-strings you created in T9 adjusted to the number from Q4 to create the below data
Click here for hint
In the console, write?separate
and think about how you used it earlier. Perhaps you can not only specify a vector to separate into
, but also specify a function, which returns a vector?
|>
peptide_data sample_n(size = 10)
# A tibble: 10 × 18
Experiment CDR3b V_gene J_gene peptide_1 peptide_2 peptide_3 peptide_4
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 eEE243 CASSYGSLYNE… TCRBV… TCRBJ… MWSFNPET… SFNPETNIL SMWSFNPET <NA>
2 eXL31 CASSLEFSVTG… TCRBV… TCRBJ… LSPRWYFYY SPRWYFYYL <NA> <NA>
3 eHH175 CASSQLAGGPE… TCRBV… TCRBJ… AFLLFLVLI FLAFLLFLV FYLCFLAFL FYLCFLAF…
4 eEE240 CASSTDRLYGY… TCRBV… TCRBJ… KPLEFGAT… <NA> <NA> <NA>
5 eOX52 CASRPGQGGQP… TCRBV… TCRBJ… MWSFNPET… SFNPETNIL SMWSFNPET <NA>
6 eMR17 unproductive TCRBV… TCRBJ… HTTDPSFL… <NA> <NA> <NA>
7 eEE224 CASSSGYQETQ… TCRBV… TCRBJ… ELYSPIFLI LYSPIFLIV QELYSPIFL VQELYSPIF
8 eEE228 CASRILAGGRT… TCRBV… TCRBJ… KLPDDFTG… <NA> <NA> <NA>
9 eEE224 CSVGLAGGADT… TCRBV… TCRBJ… DFLEYHDVR EDFLEYHD… LEYHDVRVV LEYHDVRV…
10 eXL31 CASSSGLAGVG… TCRBV… TCRBJ… YFPLQSYGF <NA> <NA> <NA>
# ℹ 10 more variables: peptide_5 <chr>, peptide_6 <chr>, peptide_7 <chr>,
# peptide_8 <chr>, peptide_9 <chr>, peptide_10 <chr>, peptide_11 <chr>,
# peptide_12 <chr>, peptide_13 <chr>, n_peptides <dbl>
Q5: Now, presumable you got a warning, discuss in your group why that is?
Q6: With respect to
peptide_n
, discuss in your group, if this is wide- or long-data?
Now, finally we will use the what we prepared for today, data pivoting. There are two functions, namely pivot_wider()
and pivot_longer()
. Also, now, we will use a trick when developing ones data pipeline, while working with new functions, that on might not be completely comfortable with. You have seen the sample_n()
function several times above and we can use that to randomly sample n
observations from data. This we can utilise to work with a smaller data set in the development face and once we are ready, we can increase this n
gradually to see if everything continues to work as anticipated.
T11: Using the
peptide_data
, run a fewsample_n()
calls with varying degree ofn
to make sure, that you get a feeling for what is going onT12: From the
peptide_data
data above, with peptide_1, peptide_2, etc. create this data set using one of the data pivoting functions. Remember to start initially with sampling a smaller data set and then work on that first! Also, once you’re sure you’re good to go, reuse thepeptide_data
variable as we don’t want huge redundant data sets floating around in our environment
Click here for hint
If the pivoting is not clear at all, then do what I do, create some example data:
<- tibble(
my_data id = str_c("id_", 1:10),
var_1 = round(rnorm(10),1),
var_2 = round(rnorm(10),1),
var_3 = round(rnorm(10),1))
…and then play around with that. A small set like the one above is easy to handle, so perhaps start with that and then pivot back and forth a few times using pivot_wider()
/pivot_longer()
. Use View()
to inspect and get a better overview of the results of pivoting.
|>
peptide_data sample_n(10)
# A tibble: 10 × 7
Experiment CDR3b V_gene J_gene n_peptides peptide_n peptide
<chr> <chr> <chr> <chr> <dbl> <chr> <chr>
1 eXL27 CASSAGTTYEQYF TCRBV05-06 TCRBJ0… 2 peptide_3 <NA>
2 eXL30 CASSSLAGGLDDEQFF TCRBV06-X TCRBJ0… 10 peptide_1 FKVSIW…
3 eHH175 CASSQGARSTEAFF TCRBV06-05 TCRBJ0… 11 peptide_3 FYLCFL…
4 eEE240 CASSEDGGNTEAFF TCRBV06-01 TCRBJ0… 4 peptide_7 <NA>
5 eOX54 CAISAPPRGGPNTGELFF TCRBV10-03 TCRBJ0… 11 peptide_6 IELSLI…
6 eMR13 CASRETPGFETQYF TCRBV06-05 TCRBJ0… 1 peptide_1 HTTDPS…
7 ePD83 CASSIATGLAGEQYF TCRBV19-01 TCRBJ0… 3 peptide_4 <NA>
8 eAV91 CASSLGLAGEAGANVLTF TCRBV07-06 TCRBJ0… 1 peptide_1 FVDGVP…
9 eEE226 CASSEAGGQSDTQYF TCRBV07-03 TCRBJ0… 7 peptide_2 GYINVF…
10 eEE226 CASSPGDRAYEQYF TCRBV07-06 TCRBJ0… 1 peptide_… <NA>
Q7: You will see some
NA
s in thepeptide
variable, discuss in your group from where these arise?Q8: How many rows and columns now and how does this compare with Q3? Discuss why/why not it is different?
T13: Now, lose the redundant variables
n_peptides
andpeptide_n
, get rid of theNA
s in thepeptide
column, and make sure that we only have unique observations (i.e. there are no repeated rows/observations).
|>
peptide_data sample_n(10)
# A tibble: 10 × 5
Experiment CDR3b V_gene J_gene peptide
<chr> <chr> <chr> <chr> <chr>
1 eOX43 CASSRRGSYEQYF TCRBV13-01 TCRBJ02-07 CFTNVYADSF
2 eQD121 CASSPRDSGLTGELFF TCRBV27-01 TCRBJ02-02 HTTDPSFLGRY
3 eOX52 CASSEGTYEQYF TCRBV02-01 TCRBJ02-07 ALALLLLDR
4 eXL30 CASSPLTSGYNEQFF TCRBV13-01 TCRBJ02-01 LLFLVLIML
5 eEE240 CASGLEPGELFF TCRBV12-05 TCRBJ02-02 FLNGSCGSV
6 eEE224 CASTEDGMNTEAFF TCRBV11-02 TCRBJ01-01 KLSYGIATV
7 eOX52 CASSLTGVYGELFF TCRBV12-X TCRBJ02-02 NVFAFPFTI
8 eOX49 CATSTADGREGSTDTQYF TCRBV15-01 TCRBJ02-03 YLCFLAFLL
9 eQD123 CASSDRTGTVNTEAFF TCRBV09-01 TCRBJ01-01 AYKTFPPTEPK
10 eEE226 CASSPWDRAYEQYF TCRBV09-01 TCRBJ02-07 KLSYGIATV
- Q8: Now how many rows and columns and is this data tidy? Discuss in your group why/why not?
Again, we turn to the stringr
package, as we need to make sure that the sequence data does indeed only contain valid characters. There are a total of 20 proteogenic amino acids, which we symbolise using ARNDCQEGHILKMFPSTWYV
.
- T14: Use the
str_detect()
function tofilter
theCDR3b
andpeptide
variables using apattern
of[^ARNDCQEGHILKMFPSTWYV]
and then play with thenegate
parameter so see what happens
Click here for hint
Again, try to play a bit around with the function in the console, type e.g.str_detect(string = "ARND", pattern = "A")
and str_detect(string = "ARND", pattern = "C")
and then recall, that the filter()
function requires a logical vector, i.e. a vector of TRUE
and FALSE
to filter the rows
- T15: Add two new variables to the data,
k_CDR3b
andk_peptide
each signifying the length of the respective sequences
Click here for hint
Again, we’re working with strings, so perhaps there is a package of interest and perhaps in that package, there is a function, which can get the length of a string?|>
peptide_data sample_n(10)
# A tibble: 10 × 7
Experiment CDR3b V_gene J_gene peptide k_CDR3b k_peptide
<chr> <chr> <chr> <chr> <chr> <int> <int>
1 eEE240 CASSLGAGSGVFF TCRBV07-02 TCRBJ01-04 VYFLQS… 13 10
2 eEE224 CSVETDTQYF TCRBV29-01 TCRBJ02-03 YLQPRT… 10 9
3 eOX52 CASSASGNPDEQFF TCRBV09-01 TCRBJ02-01 FLWLLW… 14 9
4 eAV93 CASSPGVARGGETQYF TCRBV11-03 TCRBJ02-05 LWPVTL… 16 9
5 eEE240 CASSASGQVSDTQYF TCRBV09-01 TCRBJ02-03 FVDGVP… 15 9
6 eHO141 CASSQREELITGELFF TCRBV04-01 TCRBJ02-02 KTFPPT… 16 9
7 eHO136 CASSYSYSFSSERDTQYF TCRBV06-X TCRBJ02-03 FLAFLL… 18 9
8 eHO138 CASSPEPRRRETQYF TCRBV11-02 TCRBJ02-05 MVMCGG… 15 10
9 eEE226 CASSSTSNQPQHF TCRBV11-02 TCRBJ01-05 DGVYFA… 13 10
10 eOX52 CASSQVGEVRNEQFF TCRBV04-03 TCRBJ02-01 IELSLI… 15 10
- T16: Re-create this plot
Q9: What is the most predominant length of the CDR3b-sequences?
T17: Re-create this plot
Q10: What is the most predominant length of the peptide-sequences?
Q11: Discuss in your group, if this data set is tidy or not?
|>
peptide_data sample_n(10)
# A tibble: 10 × 7
Experiment CDR3b V_gene J_gene peptide k_CDR3b k_peptide
<chr> <chr> <chr> <chr> <chr> <int> <int>
1 eXL31 CASWEGQETQYF TCRBV02-01 TCRBJ02-05 LIDFYLC… 12 9
2 eHO126 CASSLAGSEQFF TCRBV05-01 TCRBJ02-01 LSPRWYF… 12 9
3 eOX54 CSARDFKGRSYEQYF TCRBV20-X TCRBJ02-07 TLATCEL… 15 10
4 eQD126 CASSLESAMNTEAFF TCRBV07-09 TCRBJ01-01 SPRWYFY… 15 9
5 eJL158 CASTSSVNEQFF TCRBV27-01 TCRBJ02-01 SIIAYTM… 12 9
6 eOX54 CATQTLRDSLSYNEQFF TCRBV06-X TCRBJ02-01 YEQYIKW… 17 10
7 eOX43 CSADDASSGRRNTQYF TCRBV20-X TCRBJ02-03 LLFLVLI… 16 9
8 eQD108 CASSLGRTNTGELFF TCRBV07-08 TCRBJ02-02 LYSPIFL… 15 9
9 eOX52 CASSQDYEGAPPHEQYF TCRBV04-01 TCRBJ02-07 GYINVFA… 17 10
10 eEE224 CASSLGTGINEQFF TCRBV27-01 TCRBJ02-01 MEVTPSG… 14 10
Creating one data set from two data sets
Before we move onto using the family of *_join()
functions you prepared for today, we will just take a quick peek at the meta data again:
|>
meta_data sample_n(10)
# A tibble: 10 × 11
Experiment Cohort Age Gender Race A1 A2 B1 B2 C1 C2
<chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 eAV88 Healthy (N… 24 M White "A*0… "A*0… "B*2… "B*4… "C*0… "C*0…
2 eOX43 Healthy (N… 24 M White "A*0… "A*0… "B*2… "B*4… "C*0… "C*0…
3 eEE224 Healthy (N… 24 M White "A*0… "A*0… "B*2… "B*4… "C*0… "C*0…
4 eLH51 COVID-19-C… 55 M Asian "A*2… "A*3… "B*1… "B*5… "C*0… "C*1…
5 eMR26 COVID-19-C… 62 M <NA> "A*0… "A*0… "B*0… "B*5… "C*0… "C*0…
6 ePD86 COVID-19-C… 58 M White "A*0… "A*2… "B*4… "B*5… "C*0… "C*1…
7 eHO141 COVID-19-A… NA <NA> <NA> "" "" "" "" "" ""
8 eQD128 COVID-19-C… 53 F Asian "A*0… "A*1… "B*3… "B*4… "C*0… "C*0…
9 eJL164 COVID-19-B… 33 M White "A*0… "A*2… "B*1… "B*4… "C*0… "C*1…
10 eQD125 COVID-19-C… 44 M <NA> "A*0… "A*1… "B*1… "B*5… "C*0… "C*0…
Remember you can scroll in the data.
- Q12: Discuss in your group, if this data with respect to the
A1
,A2
,B1
,B2
,C1
andC2
variables is a wide or a long data format?
As with the peptide_data
, we will now have to use data pivoting again. I.e.:
- T18: use either
pivot_wider()
orpivot_longer()
to create the following data:
|>
meta_data sample_n(10)
# A tibble: 10 × 7
Experiment Cohort Age Gender Race Gene Allele
<chr> <chr> <dbl> <chr> <chr> <chr> <chr>
1 eQD112 COVID-19-Convalescent 65 M <NA> A1 A*24:02:01
2 eLH44 COVID-19-Convalescent 61 F <NA> B1 B*14:02:01
3 ePD85 Healthy (No known exposure) 27 F <NA> A2 A*29:01
4 eEE217 Healthy (No known exposure) 32 F White C2 C*05:01
5 eOX49 Healthy (No known exposure) 21 M White C1 C*03:03
6 eXL31 Healthy (No known exposure) 28 M White C2 C*16:01
7 eXL32 Healthy (No known exposure) 37 F White A1 A*01:01
8 eEE243 Healthy (No known exposure) 32 F <NA> A1 A*02:01
9 ePD82 COVID-19-Convalescent 60 F <NA> C1 C*08:01:01
10 eQD131 COVID-19-Exposed NA <NA> <NA> B1 B*15:01:01
Remember, what we are aiming for here, is to create one data set from two. So:
- Q13: Discuss in your group, which variable(s?) define the same observations between the
peptide_data
and themeta_data
?
Once you have agreed upon Experiment
, then use that knowledge to subset the meta_data
to the variables-of-interest:
|>
meta_data sample_n(10)
# A tibble: 10 × 2
Experiment Allele
<chr> <chr>
1 eLH43 "B*44:02:01"
2 eJL157 "C*07:02:01"
3 eDH113 "B*44:03"
4 eHO138 ""
5 eQD116 "B*35:01:01"
6 eMR25 ""
7 ePD87 "B*07:02:01"
8 eQD108 "A*68:01:02"
9 eHH173 "B*35:01"
10 eMR21 "A*29:02:01"
Use the View()
function again, to look at the meta_data
. Notice something? Some alleles are e.g. A*11:01
, whereas others are B*51:01:02
. You can find information on why, by visiting Nomenclature for Factors of the HLA System.
Long story short, we only want to include Field 1
(allele group) and Field 2
(Specific HLA protein). You have prepared the stringr
package for today. See if you can find a way to reduce e.g. B*51:01:02
to B*51:01
and then create a new variable Allele_F_1_2
accordingly, while also removing the ...x
(where x
is a number) subscripts from the Gene
variable (It is an artifact from having the data in a wide format, where you cannot have two variables with the same name) and also, remove any NA
s and ""
s, denoting empty entries.
Click here for hint
There are several ways this can be achieved, the easiest being to consider if perhaps a part of the string based on indices could be of interest. This term “a part of a string” is called a substring, perhaps thestringr
package contains a function work with substring? In the console, type stringr::
and hit tab
. This will display the functions available in the stringr
package. Scroll down and find the functionst starting with str_
and look for on, which might be relevant and remember you can use ?function_name
to get more information on how a given function works.
- T19: Create the following data, according to specifications above:
|>
meta_data sample_n(10)
# A tibble: 10 × 3
Experiment Allele Allele_F_1_2
<chr> <chr> <chr>
1 eQD115 C*07:02:01 C*07:02
2 eHH173 B*44:02 B*44:02
3 eQD114 A*24:02:01 A*24:02
4 eOX52 C*04:82 C*04:82
5 ePD80 A*02:01:01 A*02:01
6 eQD134 C*04:03:01 C*04:03
7 eJL154 A*29:02:01 A*29:02
8 eHO135 A*31:01:02 A*31:01
9 eMR13 C*03:04:01 C*03:04
10 eQD124 B*15:01:01 B*15:01
The asterisk, i.e. *
is a rather annoying character because of ambiguity, so:
- T20: Clean the data a bit more, by removing the asterisk and redundant variables:
|>
meta_data sample_n(size = 10)
# A tibble: 10 × 2
Experiment Allele
<chr> <chr>
1 eHH173 C05:01
2 eJL152 A24:02
3 eJL161 A02:01
4 eHO126 B57:01
5 eGK111 C07:01
6 eLH58 B52:01
7 eQD112 B35:02
8 eQD136 C03:03
9 eHO133 A32:01
10 eJL154 A29:02
Click here for hint 1
Again, thestringr
package may come in handy. Perhaps there is a function remove
, one or more such pesky characters?
Click here for hint 2
Getting a weird error? Recall, that character ambiguity needs to be “escaped”, you did this somehow earlier on…Recall the peptide_data
?
|>
peptide_data sample_n(10)
# A tibble: 10 × 7
Experiment CDR3b V_gene J_gene peptide k_CDR3b k_peptide
<chr> <chr> <chr> <chr> <chr> <int> <int>
1 eEE228 CASSLPENSYEQYF TCRBV07-08 TCRBJ02… YLDAYN… 14 9
2 eAV93 CASSLIGTSGGYNEQFF TCRBV28-01 TCRBJ02… VQELYS… 17 9
3 eXL30 CSADTGDQTQYF TCRBV20-X TCRBJ02… YLCFLA… 12 9
4 eXL30 CASSPRDITYEQYF TCRBV11-01 TCRBJ02… VLPFND… 14 10
5 eOX54 CASGTTGGAETQYF TCRBV28-01 TCRBJ02… LEYHDV… 14 9
6 eXL30 CASSLSPGLAGGWGTDTQYF TCRBV12-X TCRBJ02… LLFLVL… 20 9
7 eOX46 CASSAGPRSPYEQYF TCRBV13-01 TCRBJ02… FYLCFL… 15 10
8 eEE228 CATSRGLGTGRTDTQYF TCRBV15-01 TCRBJ02… LIDFYL… 17 9
9 eEE240 CSATSMNTEAFF TCRBV20-X TCRBJ01… FYLCFL… 12 9
10 eEE224 CSVETPSGGIHEQFF TCRBV29-01 TCRBJ02… VQELYS… 15 10
- T21: Create a
dplyr
pipeline, starting with thepeptide_data
, which joins it with themeta_data
and remember to make sure that you get only unqiue observations of rows. Save this data into a new variable namespeptide_meta_data
(If you get a warning, discuss in your group what it means?)
Click here for hint 1
Which family of functions do we use to join data? Also, perhaps here it would be prudent to start with working on a smaller data set, recall we could sample a number of rows yielding a smaller development data set
Click here for hint 2
You should get a data set of around +3.000.000, take a moment to consider how that would have been to work with in Excel? Also, in case the servers are not liking this, you can consider subsetting thepeptide_data
prior to joining to e.g. 100,000 or 10,000 rows.
|>
peptide_meta_data sample_n(10)
# A tibble: 10 × 8
Experiment CDR3b V_gene J_gene peptide k_CDR3b k_peptide Allele
<chr> <chr> <chr> <chr> <chr> <int> <int> <chr>
1 eXL30 CSVEVNSYEQYF TCRBV29… TCRBJ… AFLLFL… 12 9 C07:02
2 eOX49 CASTLWVHPDTQYF TCRBV07… TCRBJ… FYLCFL… 14 9 C03:03
3 eEE226 CASSQDPGQGVLYEQYF TCRBV04… TCRBJ… IELSLI… 17 10 B35:02
4 eAV93 CASSQEPAAGLNSPLHF TCRBV04… TCRBJ… GYINVF… 17 10 C04:01
5 eXL37 CASSGTGRSGANVLTF TCRBV07… TCRBJ… VLPFND… 16 10 C03:04
6 eHH175 CASSLDSGWSAQYF TCRBV11… TCRBJ… LLFLVL… 14 9 C16:01
7 eXL31 CASSLSSVSMAFF TCRBV27… TCRBJ… QLMCQP… 13 10 B44:03
8 eEE240 CASSALVGGYNEQFF TCRBV09… TCRBJ… FLNGSC… 15 9 C03:04
9 eEE240 CSGSGSSYEQYF TCRBV29… TCRBJ… AFLLFL… 12 9 C06:02
10 eXL30 CASSPGDTQYF TCRBV27… TCRBJ… TVYSHL… 11 9 A02:01
Analysis
Now, that we have the data in a prepared and ready-to-analyse format, let us return to the two burning questions we had:
- What characterises the peptides binding to the HLAs?
- What characterises T-cell Receptors binding to the pMHC-complexes?
Peptides binding to HLA
As we have touched upon multiple times, R
is very flexible and naturally you can also create sequence logos. Finally, let us create a binding motif using the package ggseqlogo
(More info here).
- T22: Subset the final
peptide_meta_data
data toA02:01
and unique observations of peptides of length 9 and re-create the below sequence logo
Click here for hint
You can pipe a vector of peptides intoggseqlogo
, but perhaps you first need to pull
that vector from the relevant variable in your tibble? Also, consider before that, that you’ll need to make sure, you are only looking at peptides of length 9
- T23: Repeat for e.g.
B07:02
or another of your favourite alleles
Now, let’s take a closer look at the sequence logo:
- Q14: Which positions in the peptide determines binding to HLA?
Click here for hint
Recall your Introduction to Bioinformatics course? And/or perhaps ask your fellow group members if they know?CDR3b-sequences binding to pMHC
- T24: Subset the
peptide_meta_data
, such that the length of the CDR3b is 15, the allele is A02:01 and the peptide is LLFLVLIML and re-create the below sequence logo of the CDR3b sequences:
Q15: In your group, discuss what you see?
T25: Play around with other combinations of
k_CDR3b
,Allele
, andpeptide
and inspect how the logo changes
Disclaimer: In this data set, we only get: A given CDR3b was found to recognise a given peptide in a given subject and that subject had a given haplotype - Something’s missing… Perhaps if you have had immunology, then you can spot it? There is a trick to get around this missing information, but that’s beyond scope of what we’re working with here.
Epilogue
That’s it for today - I know this is overwhelming now, but commit to it and you WILL be plenty rewarded! I hope today was at least a glimpse into the flexibility and capabilities of using tidyverse
for applied Bio Data Science
…also, noticed something? We spend maybe 80% of the time here on dealing with data-wrangling and then once we’re good to go, the analysis wasn’t that time consuming - That’s often the way it ends up going. You’ll spend a lot of time on data handling, and getting the tidyverse toolbox in your tool belt will allow you to be so much more efficient in your data wrangling, so you can get to the fun part as quickly as possible!
Today’s Assignment
After today, we are halfway through the labs of the course, so now is a good time to spend some time recalling what we have been over and practising writing a reproducible Quarto-report.
Your group assignment today is to condense the exercises into a group micro-report! Talk together and figure out how to distil the exercises from today into one small end-to-end runnable reproducible micro-report. DO NOT include ALL of the exercises, but rather include as few steps as possible to arrive at your results. Be very concise!
But WHY? WHY are you not specifying exactly what we need to hand in? Because we are training taking independent decisions, which is crucial in applied bio data science, so take a look at the combined group code, select relevant sections and condense - If you don’t make it all the way through the exercises, then condense and present what you were able to arrive at! What do you think is central/important/indispensable? Also, these hand ins are NOT for us to evaluate you, but for you to train creating products and the get feedback on your progress!
IMPORTANT: Remember to check the ASSIGNMENT GUIDELINES
…and as always - Have fun!