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.
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)"
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()
scratch_dir <- file.path("..", "..", "scratch")
output_dir <- file.path(scratch_dir, "cytoband_status", "segments")
if(!dir.exists(output_dir)) {
dir.create(output_dir, recursive = TRUE)
}
if ( params$base_run ==0 ){
histologies_file <- file.path("..", "..", "data", "pbta-histologies.tsv")
consensus_seg_file <- file.path("..", "..", "data",
"pbta-cnv-consensus.seg.gz")
} else {
histologies_file <- file.path("..", "..", "data", "pbta-histologies-base.tsv")
consensus_seg_file <- file.path("..", "copy_number_consensus_call", "results",
"pbta-cnv-consensus.seg.gz")
}
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()
))
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.numcolumn 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 fromcnv_consensus.tsvare supplementented with “neutral” segments where no call was made. These include all non-variant regions present inref/cnv_callable.bedThe 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())
add_autosomes_df <- add_ploidy_df %>%
# x and y chromosomes must be handled differently
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"
))
# 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)
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))
output_file <- file.path("..", "..", "scratch", "consensus_seg_with_status.tsv")
write_tsv(add_status_df, output_file)
bed_status_df <- add_status_df %>%
select(chrom, loc.start, loc.end, everything()) %>%
group_by(Kids_First_Biospecimen_ID) %>%
arrange(chrom, loc.start, loc.end)
losses_bed_status_df <- bed_status_df %>%
filter(status == "loss")
gains_bed_status_df <- bed_status_df %>%
filter(status == "gain")
# make lists of data frames by sample
bed_status_list <- bed_status_df %>%
group_split()
names(bed_status_list) <- bed_status_df %>%
group_keys() %>%
pull(Kids_First_Biospecimen_ID)
bed_loss_list <- losses_bed_status_df %>%
group_split()
names(bed_loss_list) <- losses_bed_status_df %>%
group_keys() %>%
pull(Kids_First_Biospecimen_ID)
bed_gain_list <- gains_bed_status_df %>%
group_split()
names(bed_gain_list) <- gains_bed_status_df %>%
group_keys() %>%
pull(Kids_First_Biospecimen_ID)
temp <- purrr::imap(bed_status_list,
~ write_tsv(.x,
file.path(
output_dir, paste0("consensus_callable.", .y, ".bed")
),
col_names = FALSE))
temp_loss <- purrr::imap(bed_loss_list,
~ write_tsv(.x,
file.path(
output_dir, paste0("consensus_loss.", .y, ".bed")
),
col_names = FALSE))
temp_gain <- purrr::imap(bed_gain_list,
~ write_tsv(.x,
file.path(
output_dir, paste0("consensus_gain.", .y, ".bed")
),
col_names = FALSE))
sessionInfo()
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)
Matrix products: default
BLAS/LAPACK: /usr/lib/libopenblasp-r0.2.19.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2
[5] readr_1.3.1 tidyr_0.8.3 tibble_2.1.3 ggplot2_3.2.0
[9] tidyverse_1.2.1
loaded via a namespace (and not attached):
[1] Rcpp_1.0.1 cellranger_1.1.0 pillar_1.4.2 compiler_3.6.0
[5] base64enc_0.1-3 tools_3.6.0 digest_0.6.20 lubridate_1.7.4
[9] jsonlite_1.6 evaluate_0.14 nlme_3.1-140 gtable_0.3.0
[13] lattice_0.20-38 pkgconfig_2.0.2 rlang_0.4.0 cli_1.1.0
[17] rstudioapi_0.10 yaml_2.2.0 haven_2.1.1 xfun_0.8
[21] withr_2.1.2 xml2_1.2.0 httr_1.4.0 knitr_1.23
[25] generics_0.0.2 hms_0.4.2 grid_3.6.0 tidyselect_0.2.5
[29] glue_1.3.1 R6_2.4.0 readxl_1.3.1 rmarkdown_1.13
[33] modelr_0.1.4 magrittr_1.5 backports_1.1.4 scales_1.0.0
[37] htmltools_0.3.6 rvest_0.3.4 assertthat_0.2.1 colorspace_1.4-1
[41] labeling_0.3 stringi_1.4.3 lazyeval_0.2.2 munsell_0.5.0
[45] broom_0.5.2 crayon_1.3.4