The pbta-histologies.tsv file contains a tumor_ploidy column, which is tumor ploidy as inferred by ControlFreeC. The copy number information should be interpreted in the light of this information (see: current version of Data Formats file).

This notebook adds ploidy information to the consensus SEG file and adds a status column that defines gain and loss broadly. Note that the consensus SEG file does not have copy number information for sex chromosomes.

Usage

This notebook is intended to be run from the command line with the following (assumes you are in the root directory of the repository):

Rscript -e "rmarkdown::render('analyses/focal-cn-file-preparation/02-add-ploidy-consensus.Rmd', clean = TRUE)"

Set up

Libraries and functions

library(tidyverse)
── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.0     ✔ purrr   0.3.2
✔ tibble  2.1.3     ✔ dplyr   0.8.3
✔ tidyr   0.8.3     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.4.0
── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

Files and directories

scratch_dir <- file.path("..", "..", "scratch")
# TODO: the consensus SEG file is not currently in the data download -- when it
# gets included we will have to change the file path here
consensus_seg_file <- file.path("..", "copy_number_consensus_call", "results", 
                                "pbta-cnv-consensus.seg.gz")
histologies_file <- file.path("..", "..", "data", "pbta-histologies.tsv")

consensus_seg_df <- read_tsv(consensus_seg_file)
Parsed with column specification:
cols(
  ID = col_character(),
  chrom = col_character(),
  loc.start = col_double(),
  loc.end = col_double(),
  num.mark = col_logical(),
  seg.mean = col_double(),
  copy.num = col_double()
)
histologies_df <- read_tsv(histologies_file,
                           col_types = cols(
                             molecular_subtype = col_character()
                           ))

Add ploidy and status

As noted above, the sex chromosomes do not have copy number information included in the consensus SEG file, but other copy number values can be missing.

From the consensus SEG file methods:

The copy.num column is the weighted median of CNVkit segment values where they exist, or Control-FREEC values in the absence of CNVkit data. Because some software (notably GISTIC) requires all samples to have the same regions called, the copy number variants from cnv_consensus.tsv are supplementented with “neutral” segments where no call was made. These include all non-variant regions present in ref/cnv_callable.bed The neutral regions are assigned copy.num 2, except on chrX and chrY, where the copy number is left NA.

For segments with missing copy.num, what chromosomes are they on?

consensus_seg_df %>%
  filter(is.na(copy.num)) %>%
  group_by(chrom) %>%
  tally()

The majority are on the sex chromosomes, which is expected. We can remove these, as we can not add status using the ploidy information without the copy.num information. We need to add in the tumor_ploidy column from the histologies file.

add_ploidy_df <- consensus_seg_df %>%
  filter(!is.na(copy.num)) %>%
  inner_join(select(histologies_df, 
                    Kids_First_Biospecimen_ID, 
                    tumor_ploidy,
                    germline_sex_estimate), 
             by = c("ID" = "Kids_First_Biospecimen_ID")) %>%
  select(-tumor_ploidy, -germline_sex_estimate, everything())

Handle autosomes first

add_autosomes_df <- add_ploidy_df %>%
  # x and y chromosomes must be handled differently - leave in for safety
  filter(!(chrom %in% c("chrX", "chrY"))) %>%
  mutate(status = case_when(
    # when the copy number is less than inferred ploidy, mark this as a loss
    copy.num < tumor_ploidy ~ "loss",
    # if copy number is higher than ploidy, mark as a gain
    copy.num > tumor_ploidy ~ "gain",
    copy.num == tumor_ploidy ~ "neutral"
  ))

Handle sex chromosomes

# this logic is consistent with what is observed in the controlfreec file
# specifically, in samples where germline sex estimate = Female, X chromosomes
# appear to be treated the same way as autosomes
add_xy_df <- add_ploidy_df %>%
  filter(chrom %in% c("chrX", "chrY")) %>%  
  mutate(status = case_when(
    germline_sex_estimate == "Male" & copy.num < (0.5 * tumor_ploidy) ~ "loss",
    germline_sex_estimate == "Male" & copy.num > (0.5 * tumor_ploidy) ~ "gain",
    germline_sex_estimate == "Male" & copy.num == (0.5 * tumor_ploidy) ~ "neutral",
    # there are some instances where there are chromosome Y segments are in
    # samples where the germline_sex_estimate is Female
    chrom == "chrY" & germline_sex_estimate == "Female" & copy.num > 0 ~ "unknown",
    chrom == "chrY" & germline_sex_estimate == "Female" & copy.num == 0 ~ "neutral",
    # mirroring controlfreec, X treated same as autosomes
    chrom == "chrX" & germline_sex_estimate == "Female" & copy.num < tumor_ploidy ~ "loss",
    chrom == "chrX" & germline_sex_estimate == "Female" & copy.num > tumor_ploidy ~ "gain",
    chrom == "chrX" & germline_sex_estimate == "Female" & copy.num == tumor_ploidy ~ "neutral"
  ))

add_status_df <- bind_rows(add_autosomes_df, add_xy_df) %>%
  select(-germline_sex_estimate) %>%
  rename(Kids_First_Biospecimen_ID = ID)

Does seg.mean agree with status?

add_status_df %>%
  filter(!is.na(seg.mean)) %>%
  ggplot(aes(x = status, y = seg.mean, group = status)) +
  geom_jitter()

add_status_df %>%
  filter(!is.na(seg.mean)) %>%
  group_by(status) %>%
  summarize(mean = mean(seg.mean), sd = sd(seg.mean))

Write to file

output_file <- file.path("..", "..", "scratch", "consensus_seg_with_status.tsv")
write_tsv(add_status_df, output_file)
LS0tCnRpdGxlOiAiQWRkIHBsb2lkeSBjb2x1bW4sIHN0YXR1cyB0byBjb25zZW5zdXMgU0VHIGZpbGUiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRvYzogVFJVRQogICAgdG9jX2Zsb2F0OiBUUlVFCmF1dGhvcjogQ2hhbnRlIEJldGhlbGwgYW5kIEphY2x5biBUYXJvbmkgZm9yIEFMU0YgQ0NETApkYXRlOiAyMDIwCi0tLQoKVGhlIGBwYnRhLWhpc3RvbG9naWVzLnRzdmAgZmlsZSBjb250YWlucyBhIGB0dW1vcl9wbG9pZHlgIGNvbHVtbiwgd2hpY2ggaXMgdHVtb3IgcGxvaWR5IGFzIGluZmVycmVkIGJ5IENvbnRyb2xGcmVlQy4KVGhlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIHNob3VsZCBiZSBpbnRlcnByZXRlZCBpbiB0aGUgbGlnaHQgb2YgdGhpcyBpbmZvcm1hdGlvbiAoc2VlOiBbY3VycmVudCB2ZXJzaW9uIG9mIERhdGEgRm9ybWF0cyBmaWxlXShodHRwczovL2dpdGh1Yi5jb20vQWxleHNMZW1vbmFkZS9PcGVuUEJUQS1hbmFseXNpcy9ibG9iLzBlNjQyZWYwNjAyYWQwNGI4ZDBmYmQ4MGJlNjI2ZDQzNzRmZWU5MjgvZG9jL2RhdGEtZm9ybWF0cy5tZCNhLW5vdGUtb24tcGxvaWR5KSkuCgpUaGlzIG5vdGVib29rIGFkZHMgcGxvaWR5IGluZm9ybWF0aW9uIHRvIHRoZSBjb25zZW5zdXMgU0VHIGZpbGUgYW5kIGFkZHMgYSBzdGF0dXMgY29sdW1uIHRoYXQgZGVmaW5lcyBnYWluIGFuZCBsb3NzIGJyb2FkbHkuCioqTm90ZSoqIHRoYXQgdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSBkb2VzIG5vdCBoYXZlIGNvcHkgbnVtYmVyIGluZm9ybWF0aW9uIGZvciBzZXggY2hyb21vc29tZXMuCgojIyBVc2FnZQoKVGhpcyBub3RlYm9vayBpcyBpbnRlbmRlZCB0byBiZSBydW4gZnJvbSB0aGUgY29tbWFuZCBsaW5lIHdpdGggdGhlIGZvbGxvd2luZyAoYXNzdW1lcyB5b3UgYXJlIGluIHRoZSByb290IGRpcmVjdG9yeSBvZiB0aGUgcmVwb3NpdG9yeSk6CgpgYGAKUnNjcmlwdCAtZSAicm1hcmtkb3duOjpyZW5kZXIoJ2FuYWx5c2VzL2ZvY2FsLWNuLWZpbGUtcHJlcGFyYXRpb24vMDItYWRkLXBsb2lkeS1jb25zZW5zdXMuUm1kJywgY2xlYW4gPSBUUlVFKSIKYGBgCgojIyBTZXQgdXAKCiMjIyBMaWJyYXJpZXMgYW5kIGZ1bmN0aW9ucwoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpgYGAKCiMjIyBGaWxlcyBhbmQgZGlyZWN0b3JpZXMKCmBgYHtyfQpzY3JhdGNoX2RpciA8LSBmaWxlLnBhdGgoIi4uIiwgIi4uIiwgInNjcmF0Y2giKQpgYGAKCmBgYHtyfQojIFRPRE86IHRoZSBjb25zZW5zdXMgU0VHIGZpbGUgaXMgbm90IGN1cnJlbnRseSBpbiB0aGUgZGF0YSBkb3dubG9hZCAtLSB3aGVuIGl0CiMgZ2V0cyBpbmNsdWRlZCB3ZSB3aWxsIGhhdmUgdG8gY2hhbmdlIHRoZSBmaWxlIHBhdGggaGVyZQpjb25zZW5zdXNfc2VnX2ZpbGUgPC0gZmlsZS5wYXRoKCIuLiIsICJjb3B5X251bWJlcl9jb25zZW5zdXNfY2FsbCIsICJyZXN1bHRzIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgInBidGEtY252LWNvbnNlbnN1cy5zZWcuZ3oiKQpoaXN0b2xvZ2llc19maWxlIDwtIGZpbGUucGF0aCgiLi4iLCAiLi4iLCAiZGF0YSIsICJwYnRhLWhpc3RvbG9naWVzLnRzdiIpCgpjb25zZW5zdXNfc2VnX2RmIDwtIHJlYWRfdHN2KGNvbnNlbnN1c19zZWdfZmlsZSkKaGlzdG9sb2dpZXNfZGYgPC0gcmVhZF90c3YoaGlzdG9sb2dpZXNfZmlsZSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgY29sX3R5cGVzID0gY29scygKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtb2xlY3VsYXJfc3VidHlwZSA9IGNvbF9jaGFyYWN0ZXIoKQogICAgICAgICAgICAgICAgICAgICAgICAgICApKQpgYGAKCiMjIEFkZCBwbG9pZHkgYW5kIHN0YXR1cwoKQXMgbm90ZWQgYWJvdmUsIHRoZSBzZXggY2hyb21vc29tZXMgZG8gbm90IGhhdmUgY29weSBudW1iZXIgaW5mb3JtYXRpb24gaW5jbHVkZWQgaW4gdGhlIGNvbnNlbnN1cyBTRUcgZmlsZSwgYnV0IG90aGVyIGNvcHkgbnVtYmVyIHZhbHVlcyBjYW4gYmUgbWlzc2luZy4KCkZyb20gdGhlIFtjb25zZW5zdXMgU0VHIGZpbGUgbWV0aG9kc10oaHR0cHM6Ly9naXRodWIuY29tL0FsZXhzTGVtb25hZGUvT3BlblBCVEEtYW5hbHlzaXMvdHJlZS8wZTY0MmVmMDYwMmFkMDRiOGQwZmJkODBiZTYyNmQ0Mzc0ZmVlOTI4L2FuYWx5c2VzL2NvcHlfbnVtYmVyX2NvbnNlbnN1c19jYWxsI2NvbnNlbnN1cy1jbnYtY3JlYXRpb24pOiAKCj4gVGhlIGBjb3B5Lm51bWAgY29sdW1uIGlzIHRoZSB3ZWlnaHRlZCBtZWRpYW4gb2YgQ05Wa2l0IHNlZ21lbnQgdmFsdWVzIHdoZXJlIHRoZXkgZXhpc3QsIG9yIENvbnRyb2wtRlJFRUMgdmFsdWVzIGluIHRoZSBhYnNlbmNlIG9mIENOVmtpdCBkYXRhLiBCZWNhdXNlIHNvbWUgc29mdHdhcmUgKG5vdGFibHkgR0lTVElDKSByZXF1aXJlcyBhbGwgc2FtcGxlcyB0byBoYXZlIHRoZSBzYW1lIHJlZ2lvbnMgY2FsbGVkLCB0aGUgY29weSBudW1iZXIgdmFyaWFudHMgZnJvbSBgY252X2NvbnNlbnN1cy50c3ZgIGFyZSBzdXBwbGVtZW50ZW50ZWQgd2l0aCAibmV1dHJhbCIgc2VnbWVudHMgd2hlcmUgbm8gY2FsbCB3YXMgbWFkZS4gVGhlc2UgaW5jbHVkZSBhbGwgbm9uLXZhcmlhbnQgcmVnaW9ucyBwcmVzZW50IGluIGByZWYvY252X2NhbGxhYmxlLmJlZGAgVGhlIG5ldXRyYWwgcmVnaW9ucyBhcmUgYXNzaWduZWQgY29weS5udW0gMiwgZXhjZXB0IG9uIGNoclggYW5kIGNoclksIHdoZXJlIHRoZSBjb3B5IG51bWJlciBpcyBsZWZ0IE5BLgoKRm9yIHNlZ21lbnRzIHdpdGggbWlzc2luZyBgY29weS5udW1gLCB3aGF0IGNocm9tb3NvbWVzIGFyZSB0aGV5IG9uPwoKYGBge3J9CmNvbnNlbnN1c19zZWdfZGYgJT4lCiAgZmlsdGVyKGlzLm5hKGNvcHkubnVtKSkgJT4lCiAgZ3JvdXBfYnkoY2hyb20pICU+JQogIHRhbGx5KCkKYGBgCgpUaGUgbWFqb3JpdHkgYXJlIG9uIHRoZSBzZXggY2hyb21vc29tZXMsIHdoaWNoIGlzIGV4cGVjdGVkLgpXZSBjYW4gcmVtb3ZlIHRoZXNlLCBhcyB3ZSBjYW4gbm90IGFkZCBzdGF0dXMgdXNpbmcgdGhlIHBsb2lkeSBpbmZvcm1hdGlvbiB3aXRob3V0IHRoZSBgY29weS5udW1gIGluZm9ybWF0aW9uLgpXZSBuZWVkIHRvIGFkZCBpbiB0aGUgYHR1bW9yX3Bsb2lkeWAgY29sdW1uIGZyb20gdGhlIGhpc3RvbG9naWVzIGZpbGUuCgpgYGB7cn0KYWRkX3Bsb2lkeV9kZiA8LSBjb25zZW5zdXNfc2VnX2RmICU+JQogIGZpbHRlcighaXMubmEoY29weS5udW0pKSAlPiUKICBpbm5lcl9qb2luKHNlbGVjdChoaXN0b2xvZ2llc19kZiwgCiAgICAgICAgICAgICAgICAgICAgS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCwgCiAgICAgICAgICAgICAgICAgICAgdHVtb3JfcGxvaWR5LAogICAgICAgICAgICAgICAgICAgIGdlcm1saW5lX3NleF9lc3RpbWF0ZSksIAogICAgICAgICAgICAgYnkgPSBjKCJJRCIgPSAiS2lkc19GaXJzdF9CaW9zcGVjaW1lbl9JRCIpKSAlPiUKICBzZWxlY3QoLXR1bW9yX3Bsb2lkeSwgLWdlcm1saW5lX3NleF9lc3RpbWF0ZSwgZXZlcnl0aGluZygpKQpgYGAKCiMjIyBIYW5kbGUgYXV0b3NvbWVzIGZpcnN0CgpgYGB7cn0KYWRkX2F1dG9zb21lc19kZiA8LSBhZGRfcGxvaWR5X2RmICU+JQogICMgeCBhbmQgeSBjaHJvbW9zb21lcyBtdXN0IGJlIGhhbmRsZWQgZGlmZmVyZW50bHkgLSBsZWF2ZSBpbiBmb3Igc2FmZXR5CiAgZmlsdGVyKCEoY2hyb20gJWluJSBjKCJjaHJYIiwgImNoclkiKSkpICU+JQogIG11dGF0ZShzdGF0dXMgPSBjYXNlX3doZW4oCiAgICAjIHdoZW4gdGhlIGNvcHkgbnVtYmVyIGlzIGxlc3MgdGhhbiBpbmZlcnJlZCBwbG9pZHksIG1hcmsgdGhpcyBhcyBhIGxvc3MKICAgIGNvcHkubnVtIDwgdHVtb3JfcGxvaWR5IH4gImxvc3MiLAogICAgIyBpZiBjb3B5IG51bWJlciBpcyBoaWdoZXIgdGhhbiBwbG9pZHksIG1hcmsgYXMgYSBnYWluCiAgICBjb3B5Lm51bSA+IHR1bW9yX3Bsb2lkeSB+ICJnYWluIiwKICAgIGNvcHkubnVtID09IHR1bW9yX3Bsb2lkeSB+ICJuZXV0cmFsIgogICkpCmBgYAoKIyMjIEhhbmRsZSBzZXggY2hyb21vc29tZXMKCmBgYHtyfQojIHRoaXMgbG9naWMgaXMgY29uc2lzdGVudCB3aXRoIHdoYXQgaXMgb2JzZXJ2ZWQgaW4gdGhlIGNvbnRyb2xmcmVlYyBmaWxlCiMgc3BlY2lmaWNhbGx5LCBpbiBzYW1wbGVzIHdoZXJlIGdlcm1saW5lIHNleCBlc3RpbWF0ZSA9IEZlbWFsZSwgWCBjaHJvbW9zb21lcwojIGFwcGVhciB0byBiZSB0cmVhdGVkIHRoZSBzYW1lIHdheSBhcyBhdXRvc29tZXMKYWRkX3h5X2RmIDwtIGFkZF9wbG9pZHlfZGYgJT4lCiAgZmlsdGVyKGNocm9tICVpbiUgYygiY2hyWCIsICJjaHJZIikpICU+JSAgCiAgbXV0YXRlKHN0YXR1cyA9IGNhc2Vfd2hlbigKICAgIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiTWFsZSIgJiBjb3B5Lm51bSA8ICgwLjUgKiB0dW1vcl9wbG9pZHkpIH4gImxvc3MiLAogICAgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJNYWxlIiAmIGNvcHkubnVtID4gKDAuNSAqIHR1bW9yX3Bsb2lkeSkgfiAiZ2FpbiIsCiAgICBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIk1hbGUiICYgY29weS5udW0gPT0gKDAuNSAqIHR1bW9yX3Bsb2lkeSkgfiAibmV1dHJhbCIsCiAgICAjIHRoZXJlIGFyZSBzb21lIGluc3RhbmNlcyB3aGVyZSB0aGVyZSBhcmUgY2hyb21vc29tZSBZIHNlZ21lbnRzIGFyZSBpbgogICAgIyBzYW1wbGVzIHdoZXJlIHRoZSBnZXJtbGluZV9zZXhfZXN0aW1hdGUgaXMgRmVtYWxlCiAgICBjaHJvbSA9PSAiY2hyWSIgJiBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIkZlbWFsZSIgJiBjb3B5Lm51bSA+IDAgfiAidW5rbm93biIsCiAgICBjaHJvbSA9PSAiY2hyWSIgJiBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIkZlbWFsZSIgJiBjb3B5Lm51bSA9PSAwIH4gIm5ldXRyYWwiLAogICAgIyBtaXJyb3JpbmcgY29udHJvbGZyZWVjLCBYIHRyZWF0ZWQgc2FtZSBhcyBhdXRvc29tZXMKICAgIGNocm9tID09ICJjaHJYIiAmIGdlcm1saW5lX3NleF9lc3RpbWF0ZSA9PSAiRmVtYWxlIiAmIGNvcHkubnVtIDwgdHVtb3JfcGxvaWR5IH4gImxvc3MiLAogICAgY2hyb20gPT0gImNoclgiICYgZ2VybWxpbmVfc2V4X2VzdGltYXRlID09ICJGZW1hbGUiICYgY29weS5udW0gPiB0dW1vcl9wbG9pZHkgfiAiZ2FpbiIsCiAgICBjaHJvbSA9PSAiY2hyWCIgJiBnZXJtbGluZV9zZXhfZXN0aW1hdGUgPT0gIkZlbWFsZSIgJiBjb3B5Lm51bSA9PSB0dW1vcl9wbG9pZHkgfiAibmV1dHJhbCIKICApKQoKYWRkX3N0YXR1c19kZiA8LSBiaW5kX3Jvd3MoYWRkX2F1dG9zb21lc19kZiwgYWRkX3h5X2RmKSAlPiUKICBzZWxlY3QoLWdlcm1saW5lX3NleF9lc3RpbWF0ZSkgJT4lCiAgcmVuYW1lKEtpZHNfRmlyc3RfQmlvc3BlY2ltZW5fSUQgPSBJRCkKYGBgCgojIyMgRG9lcyBgc2VnLm1lYW5gIGFncmVlIHdpdGggYHN0YXR1c2A/CgpgYGB7cn0KYWRkX3N0YXR1c19kZiAlPiUKICBmaWx0ZXIoIWlzLm5hKHNlZy5tZWFuKSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0gc3RhdHVzLCB5ID0gc2VnLm1lYW4sIGdyb3VwID0gc3RhdHVzKSkgKwogIGdlb21faml0dGVyKCkKYGBgCgpgYGB7cn0KYWRkX3N0YXR1c19kZiAlPiUKICBmaWx0ZXIoIWlzLm5hKHNlZy5tZWFuKSkgJT4lCiAgZ3JvdXBfYnkoc3RhdHVzKSAlPiUKICBzdW1tYXJpemUobWVhbiA9IG1lYW4oc2VnLm1lYW4pLCBzZCA9IHNkKHNlZy5tZWFuKSkKYGBgCgojIyMgV3JpdGUgdG8gZmlsZQoKYGBge3J9Cm91dHB1dF9maWxlIDwtIGZpbGUucGF0aCgiLi4iLCAiLi4iLCAic2NyYXRjaCIsICJjb25zZW5zdXNfc2VnX3dpdGhfc3RhdHVzLnRzdiIpCndyaXRlX3RzdihhZGRfc3RhdHVzX2RmLCBvdXRwdXRfZmlsZSkKYGBgCg==