In this activity, we will explore real clinical data from breast cancer patients in The Cancer Genome Atlas (TCGA).
The goal is not only to learn R. The goal is to use R to ask scientific questions:
Later in DREAM-High, we will connect this type of clinical information to gene expression data. That is where computational biology becomes especially powerful: we can ask how molecular patterns relate to patient and tumor characteristics.
Before editing, use:
File → Save As…
Save a copy with your name in the file name.
This gives you your own workbook while preserving the original.
This activity uses:
brca_clin.csv — a simplified TCGA breast cancer
clinical datasetduct_cells.jpg — a schematic of breast duct
changesmetastasis.JPG — a schematic of metastatic spreadThe files should be located in the shared data
directory.
# This chunk sets up file path for the activity.
data_dir <- "/shared/dreamhigh/data"
Breast cancer often begins in the ducts or lobules of the breast. As cancer develops, cells can grow abnormally, fill spaces where they do not belong, and eventually invade nearby tissue.
The figure below shows a simplified progression from a normal duct to invasive ductal cancer.
Think before moving on:
What do you notice as the duct changes from normal tissue to invasive cancer?
Your answer: The invasive cancer starts irregular cell growth that eventually leaves the duct, which allows the cancerous cells to spread beyond just one duct and into other body systems. —
A major danger of cancer is metastasis: cancer cells can leave the original tumor, move through the body, and grow in distant organs.
Cancer metastasis is a multi-step process in which tumor cells spread from the primary tumor to distant organs. - (1) Tumor cells acquire invasive properties and detach from the primary tumor. - (2) Cancer cells migrate through surrounding tissue and invade nearby blood vessels. - (3) Tumor cells survive in the circulation and travel through the bloodstream. - (4) Circulating tumor cells exit blood vessels (extravasation) at distant sites. - (5) Cancer cells colonize distant organs and form metastatic tumors.
The figure also highlights important biological processes involved in metastasis, including:
Metastasis is biologically complex and difficult to predict or treat.
In clinical cancer datasets, we often look for features that tell us something about tumor severity, spread, treatment, or patient outcome.
We will use a CSV file. CSV stands for comma-separated values. It is a plain-text table format that can be opened in R, Excel, or many other programs.
brca_clin_df <- read.csv(
file.path(data_dir, "brca_clin.csv"),
stringsAsFactors = FALSE
)
The object name brca_clin_df is descriptive:
brca = breast cancerclin = clinical datadf = data frameA data frame is a table:
dim(brca_clin_df)
## [1] 1082 27
The first number is the number of rows. The second number is the number of columns.
n_patients <- nrow(brca_clin_df)
n_features <- ncol(brca_clin_df)
paste("This dataset contains", n_patients, "patients and", n_features, "clinical features.")
## [1] "This dataset contains 1082 patients and 27 clinical features."
colnames(brca_clin_df)
## [1] "bcr_patient_barcode"
## [2] "gender"
## [3] "race"
## [4] "ethnicity"
## [5] "age_at_diagnosis"
## [6] "year_of_initial_pathologic_diagnosis"
## [7] "vital_status"
## [8] "menopause_status"
## [9] "tumor_status"
## [10] "margin_status"
## [11] "days_to_last_followup"
## [12] "prior_dx"
## [13] "new_tumor_event_after_initial_treatment"
## [14] "radiation_therapy"
## [15] "histological_type"
## [16] "pathologic_T"
## [17] "pathologic_M"
## [18] "pathologic_N"
## [19] "pathologic_stage_sub"
## [20] "pathologic_stage"
## [21] "lymph_node_examined_count"
## [22] "number_of_lymphnodes_positive"
## [23] "initial_diagnosis_method"
## [24] "surgical_procedure"
## [25] "estrogen_receptor_status"
## [26] "progesterone_receptor_status"
## [27] "her2_receptor_status"
Question: Which column names look familiar? Which ones are confusing?
Your answer:Most look familiar, like the age, gender, race, age at diagnosis, etc. The confusing ones are related to the specifics of the cancer, like the receptor status and the various pathological specifics.
head(brca_clin_df)
## bcr_patient_barcode gender race ethnicity
## 1 TCGA-3C-AAAU FEMALE WHITE NOT HISPANIC OR LATINO
## 2 TCGA-3C-AALI FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 3 TCGA-3C-AALJ FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 4 TCGA-3C-AALK FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 5 TCGA-4H-AAAK FEMALE WHITE NOT HISPANIC OR LATINO
## 6 TCGA-5L-AAT0 FEMALE WHITE HISPANIC OR LATINO
## age_at_diagnosis year_of_initial_pathologic_diagnosis vital_status
## 1 55 2004 Alive
## 2 50 2003 Alive
## 3 62 2011 Alive
## 4 52 2011 Alive
## 5 50 2013 Alive
## 6 42 2010 Alive
## menopause_status
## 1 Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement)
## 2 Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)
## 3 Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)
## 4 [Unknown]
## 5 Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)
## 6 Post (prior bilateral ovariectomy OR >12 mo since LMP with no prior hysterectomy)
## tumor_status margin_status days_to_last_followup prior_dx
## 1 WITH TUMOR Negative 4047 No
## 2 TUMOR FREE Negative 4005 No
## 3 TUMOR FREE Negative 1474 No
## 4 TUMOR FREE Close 1448 No
## 5 TUMOR FREE Negative 348 No
## 6 TUMOR FREE Positive 1477 Yes
## new_tumor_event_after_initial_treatment radiation_therapy
## 1 NO NO
## 2 NO YES
## 3 NO NO
## 4 NO NO
## 5 NO NO
## 6 NO YES
## histological_type pathologic_T pathologic_M pathologic_N
## 1 Infiltrating Lobular Carcinoma TX MX NX
## 2 Infiltrating Ductal Carcinoma T2 M0 N1a
## 3 Infiltrating Ductal Carcinoma T2 M0 N1a
## 4 Infiltrating Ductal Carcinoma T1c M0 N0 (i+)
## 5 Infiltrating Lobular Carcinoma T2 M0 N2a
## 6 Infiltrating Lobular Carcinoma T2 M0 N0
## pathologic_stage_sub pathologic_stage lymph_node_examined_count
## 1 Stage X X 13
## 2 Stage IIB II 15
## 3 Stage IIB II 23
## 4 Stage IA I 2
## 5 Stage IIIA III 14
## 6 Stage IIA II 8
## number_of_lymphnodes_positive initial_diagnosis_method
## 1 4 [Not Available]
## 2 1 Core needle biopsy
## 3 1 Core needle biopsy
## 4 0 Core needle biopsy
## 5 4 Core needle biopsy
## 6 0 Incisional Biopsy
## surgical_procedure estrogen_receptor_status
## 1 Modified Radical Mastectomy Positive
## 2 Lumpectomy Positive
## 3 Modified Radical Mastectomy Positive
## 4 Simple Mastectomy Positive
## 5 Modified Radical Mastectomy Positive
## 6 Modified Radical Mastectomy Positive
## progesterone_receptor_status her2_receptor_status
## 1 Positive Negative
## 2 Positive Positive
## 3 Positive Indeterminate
## 4 Positive Positive
## 5 Positive Equivocal
## 6 Positive Negative
For a larger view, you can run this interactively:
View(brca_clin_df)
brca_clin_df[1, 1]
## [1] "TCGA-3C-AAAU"
The first number selects the row. The second number selects the column.
t(brca_clin_df[1, ])
## 1
## bcr_patient_barcode "TCGA-3C-AAAU"
## gender "FEMALE"
## race "WHITE"
## ethnicity "NOT HISPANIC OR LATINO"
## age_at_diagnosis "55"
## year_of_initial_pathologic_diagnosis "2004"
## vital_status "Alive"
## menopause_status "Pre (<6 months since LMP AND no prior bilateral ovariectomy AND not on estrogen replacement)"
## tumor_status "WITH TUMOR"
## margin_status "Negative"
## days_to_last_followup "4047"
## prior_dx "No"
## new_tumor_event_after_initial_treatment "NO"
## radiation_therapy "NO"
## histological_type "Infiltrating Lobular Carcinoma"
## pathologic_T "TX"
## pathologic_M "MX"
## pathologic_N "NX"
## pathologic_stage_sub "Stage X"
## pathologic_stage "X"
## lymph_node_examined_count "13"
## number_of_lymphnodes_positive "4"
## initial_diagnosis_method "[Not Available]"
## surgical_procedure "Modified Radical Mastectomy"
## estrogen_receptor_status "Positive"
## progesterone_receptor_status "Positive"
## her2_receptor_status "Negative"
The function t() transposes the row so it is easier to
read.
brca_clin_df[c(1, 2, 3), 1:6]
## bcr_patient_barcode gender race ethnicity
## 1 TCGA-3C-AAAU FEMALE WHITE NOT HISPANIC OR LATINO
## 2 TCGA-3C-AALI FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## 3 TCGA-3C-AALJ FEMALE BLACK OR AFRICAN AMERICAN NOT HISPANIC OR LATINO
## age_at_diagnosis year_of_initial_pathologic_diagnosis
## 1 55 2004
## 2 50 2003
## 3 62 2011
Prediction question: What do you think this code will show?
brca_clin_df[c(10, 20, 30), c("gender", "age_at_diagnosis", "vital_status")]
## gender age_at_diagnosis vital_status
## 10 FEMALE 56 Alive
## 20 FEMALE 40 Alive
## 30 FEMALE 34 Alive
Your answer:It will show the gender, age at diagnosis, and vital status for the specific patients in rows 10, 20, 30
Clinical data are not perfectly clean. Some values are missing, unknown, not evaluated, or not available.
Let us look for values that indicate missing or uncertain information.
missing_like_values <- c("[Not Available]", "[Not Evaluated]", "[Unknown]", "Indeterminate", "[Discrepancy]", "")
missing_like_count <- numeric(ncol(brca_clin_df))
for (i in 1:ncol(brca_clin_df)) {
missing_like_count[i] <- sum(brca_clin_df[[i]] %in% missing_like_values | is.na(brca_clin_df[[i]]))
}
missing_summary <- data.frame(
feature = colnames(brca_clin_df),
missing_or_uncertain_count = missing_like_count,
percent = round(100 * missing_like_count / nrow(brca_clin_df), 1)
)
missing_summary
## feature missing_or_uncertain_count percent
## 1 bcr_patient_barcode 0 0.0
## 2 gender 0 0.0
## 3 race 90 8.3
## 4 ethnicity 169 15.6
## 5 age_at_diagnosis 1 0.1
## 6 year_of_initial_pathologic_diagnosis 2 0.2
## 7 vital_status 0 0.0
## 8 menopause_status 87 8.0
## 9 tumor_status 36 3.3
## 10 margin_status 66 6.1
## 11 days_to_last_followup 103 9.5
## 12 prior_dx 1 0.1
## 13 new_tumor_event_after_initial_treatment 198 18.3
## 14 radiation_therapy 99 9.1
## 15 histological_type 1 0.1
## 16 pathologic_T 0 0.0
## 17 pathologic_M 0 0.0
## 18 pathologic_N 0 0.0
## 19 pathologic_stage_sub 5 0.5
## 20 pathologic_stage 0 0.0
## 21 lymph_node_examined_count 121 11.2
## 22 number_of_lymphnodes_positive 163 15.1
## 23 initial_diagnosis_method 93 8.6
## 24 surgical_procedure 55 5.1
## 25 estrogen_receptor_status 50 4.6
## 26 progesterone_receptor_status 53 4.9
## 27 her2_receptor_status 190 17.6
ordered_missing <- missing_summary[order(missing_summary$missing_or_uncertain_count, decreasing = TRUE), ]
barplot(
ordered_missing$missing_or_uncertain_count,
names.arg = ordered_missing$feature,
las = 2,
cex.names = 0.7,
main = "Missing or uncertain values by clinical feature",
ylab = "Number of patients"
)
Question: Which clinical features have the most missing or uncertain information?
Your answer:The receptor status and initial treatment have the highest number of patients with missing information.
This is an important lesson: real biomedical datasets are powerful, but they are also imperfect.
A categorical variable describes groups or labels.
Examples:
The base R function table() counts how many times each
value appears.
table(brca_clin_df$gender)
##
## FEMALE MALE
## 1070 12
barplot(
table(brca_clin_df$gender),
width = 0.5,
space = 0.5,
main = "Gender distribution in the TCGA breast cancer cohort",
ylab = "Number of patients",
las = 1
)
Question: What does this tell us about the cohort?
Your answer:The cohort is mostly female, and the male cases are rare and may require additional attention as the relationships drawn based on data may be heavily based on only female cases.
table(brca_clin_df$race)
##
## [Not Available] [Not Evaluated]
## 87 3
## AMERICAN INDIAN OR ALASKA NATIVE ASIAN
## 1 61
## BLACK OR AFRICAN AMERICAN WHITE
## 183 747
race_counts <- sort(table(brca_clin_df$race), decreasing = TRUE)
barplot(
race_counts,
main = "Race distribution",
ylab = "Number of patients",
las = 2,
cex.names = 0.8
)
Discussion question: Why is it important to look at who is represented in a biomedical dataset?
Your answer:Patients from different ethnicities will have different mutations, backgrounds, and risk factors, and responses to treatments. If there is only a majority of one ethnicity in the cohort, results cannot be generalized.
A numerical variable contains numbers.
Example:
summary(brca_clin_df$age_at_diagnosis)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
## 26.00 49.00 58.00 58.39 67.00 90.00 1
hist(
brca_clin_df$age_at_diagnosis,
breaks = 20,
main = "Age at diagnosis",
xlab = "Age at diagnosis",
ylab = "Number of patients"
)
Question: Around what ages are most patients diagnosed in this dataset?
Your answer:Late 50s and beyond, mostly in the mid 50s- mid 60s range.
Breast cancer tumors are often tested for receptors that can influence growth and treatment.
Important receptor features include:
table(brca_clin_df$estrogen_receptor_status)
##
## [Not Evaluated] Indeterminate Negative Positive
## 48 2 236 796
barplot(
table(brca_clin_df$estrogen_receptor_status),
main = "Estrogen receptor status",
ylab = "Number of patients",
las = 2
)
table(brca_clin_df$progesterone_receptor_status)
##
## [Not Evaluated] Indeterminate Negative Positive
## 49 4 340 689
table(brca_clin_df$her2_receptor_status)
##
## [Not Available] [Not Evaluated] Equivocal Indeterminate Negative
## 8 170 177 12 554
## Positive
## 161
A powerful part of data analysis is selecting patients who meet more than one condition.
For example, how many patients are both ER-positive and PR-positive?
ERpos_PRpos <- brca_clin_df[
brca_clin_df$estrogen_receptor_status == "Positive" &
brca_clin_df$progesterone_receptor_status == "Positive",
]
nrow(ERpos_PRpos)
## [1] 672
paste("There are", nrow(ERpos_PRpos), "patients with both ER-positive and PR-positive tumors.")
## [1] "There are 672 patients with both ER-positive and PR-positive tumors."
Triple-negative breast cancer means the tumor is:
These cancers can be more difficult to treat because they do not have these common therapeutic targets.
tnbc <- brca_clin_df[
brca_clin_df$estrogen_receptor_status == "Negative" &
brca_clin_df$progesterone_receptor_status == "Negative" &
brca_clin_df$her2_receptor_status == "Negative",
]
nrow(tnbc)
## [1] 114
tnbc_percent <- round(100 * nrow(tnbc) / nrow(brca_clin_df), 1)
paste(
"Triple-negative breast cancer cases make up",
tnbc_percent,
"% of this dataset."
)
## [1] "Triple-negative breast cancer cases make up 10.5 % of this dataset."
So far, we summarized one variable at a time. But clinical questions often involve relationships between variables.
er_vital_table <- table(
brca_clin_df$estrogen_receptor_status,
brca_clin_df$vital_status
)
er_vital_table
##
## Alive Dead
## [Not Evaluated] 40 8
## Indeterminate 0 2
## Negative 195 41
## Positive 697 99
Add row percentages:
round(100 * prop.table(er_vital_table, margin = 1), 1)
##
## Alive Dead
## [Not Evaluated] 83.3 16.7
## Indeterminate 0.0 100.0
## Negative 82.6 17.4
## Positive 87.6 12.4
This table asks:
Within each ER status group, what percent of patients are alive or dead?
A mosaic plot is a base R visualization for two categorical variables.
rownames(er_vital_table) <- c("NE","ND","Negative","Positive")
mosaicplot(
er_vital_table,
main = "Estrogen receptor status and vital status",
xlab = "Estrogen receptor status",
ylab = "Vital status",
color = TRUE,
las = 1
)
Question: What can you conclude from this plot? What should you be careful not to overclaim?
Your answer: Patients with positive estrogen receptor status mostly are alive; patients with negative estrogen receptor status also skew to the alive status, but there is a larger chunk that is dead. From this, we should not derive overly generalizing statements like: “patients with positive estrogen receptor status are more likely to live than those with negative status”.
Histology describes what the tumor tissue looks like under a microscope.
histology_counts <- sort(table(brca_clin_df$histological_type), decreasing = TRUE)
histology_counts
##
## Infiltrating Ductal Carcinoma Infiltrating Lobular Carcinoma
## 774 201
## Other, specify Mixed Histology (please specify)
## 45 29
## Mucinous Carcinoma Metaplastic Carcinoma
## 17 8
## Medullary Carcinoma [Not Available]
## 6 1
## Infiltrating Carcinoma NOS
## 1
par(mar = c(10, 4, 4, 2))
barplot(
histology_counts,
main = "Histological type",
ylab = "Number of patients",
las = 2,
cex.names = 0.75
)
Most tumors in this dataset are infiltrating ductal carcinoma or infiltrating lobular carcinoma.
N_ductal <- histology_counts["Infiltrating Ductal Carcinoma"]
percent_ductal <- round(100 * N_ductal / nrow(brca_clin_df), 1)
paste(
percent_ductal,
"% of patients have infiltrating ductal carcinoma.",
sep = ""
)
## [1] "71.5% of patients have infiltrating ductal carcinoma."
Pathologic stage (pathologic_stage) describes how far
the cancer has progressed.
# Uncomment the next line and place the correct feature after the $ sign
# See the hint above
table(brca_clin_df$pathologic_stage)
##
## [Not vailable] I II III IV
## 5 179 615 250 19
## X
## 14
How many patients have stage II breast cancer?
# Uncomment the next two lines and place the symbol for stage II between the parentheses
N_stageII <- table(brca_clin_df$pathologic_stage)["II"]
N_stageII
## II
## 615
# Uncomment all of the lines below
percent_stageII <- round(100 * N_stageII / nrow(brca_clin_df), 1)
paste( "There are", N_stageII, "patients with stage II breast cancer, representing",percent_stageII, "% of the cohort." )
## [1] "There are 615 patients with stage II breast cancer, representing 56.8 % of the cohort."
Make a bar plot of cancer stage
(pathological_stage).
# Uncomment the next line and place the correct feature after the $ sign
stage_counts <- table(brca_clin_df$pathological_stage)
# Uncomment the next six lines. Replace X with the object you want to make the barplot for.
barplot(
percent_stageII, main = "Pathologic stage", xlab = "Stage", ylab = "Number of patients" )
Cancer staging often depends partly on whether cancer has spread to lymph nodes.
The column number_of_lymphnodes_positive looks numeric,
but because it contains values such as [Not Available], R
may read it as text.
We can convert it carefully.
lymph_nodes_positive <- suppressWarnings(as.numeric(brca_clin_df$number_of_lymphnodes_positive))
summary(lymph_nodes_positive)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
## 0.000 0.000 1.000 2.382 2.000 35.000 163
The warning is suppressed because R turns non-numeric values like
[Not Available] into NA.
Now we can compare positive lymph node counts by pathologic stage.
boxplot(
lymph_nodes_positive ~ brca_clin_df$pathologic_stage,
main = "Positive lymph nodes by pathologic stage",
xlab = "Pathologic stage",
ylab = "Number of positive lymph nodes"
)
Question: What pattern do you see?
Your answer: The median amount of lymph nodes peaks in Stage 3 patients, and falls in the later stages. Initial stages 1, 2 do not exhibit significant amount of lymph nodes, although there are more outliers to this pattern in stage 2.
In this activity, you used R to:
Most importantly, you saw that real biomedical data are both powerful and messy.
That is normal. A major part of computational biology is learning how to reason carefully with complex, imperfect data.
Save your file. Then click:
Knit → Knit to HTML
Your final HTML report should include:
This is one reason R Markdown is powerful: it combines writing, code, results, and figures in one reproducible document.
In later DREAM-High activities, we will connect these clinical features to gene expression data from the same cancer type.
That will allow us to ask a deeper systems biology question:
Can molecular data reveal patterns that help us understand cancer subtype, treatment, and disease behavior?