From this comment on #182:
may advise starting with ATRT’s and SMARCB1 deletions
Here, I’m looking at the output of 01-prepare-cn-file.R to see if we recover that information in the segment to gene mapping. We’re going to get the data from a pull request/branch that is currently open. This will get whatever the current files are for that branch.
wget --quiet -N -P results https://raw.githubusercontent.com/jaclyn-taroni/OpenPBTA-analysis/186-both-formats/analyses/focal-cn-file-preparation/results/cnvkit_annotated_cn_autosomes.tsv.gz
wget --quiet -N -P results https://raw.githubusercontent.com/jaclyn-taroni/OpenPBTA-analysis/186-both-formats/analyses/focal-cn-file-preparation/results/controlfreec_annotated_cn_autosomes.tsv.gz `%>%` <- dplyr::`%>%`Registered S3 method overwritten by 'dplyr':
method from
print.rowwise_df # read in the output of `01-prepare-cn-file.R`
cnvkit_df <- readr::read_tsv(file.path("results",
"cnvkit_annotated_cn_autosomes.tsv.gz"))Parsed with column specification:
cols(
biospecimen_id = [31mcol_character()[39m,
status = [31mcol_character()[39m,
copy_number = [32mcol_double()[39m,
ploidy = [32mcol_double()[39m,
gene_symbol = [31mcol_character()[39m,
cytoband = [31mcol_character()[39m,
ensembl = [31mcol_character()[39m
)controlfreec_df <- readr::read_tsv(file.path("results",
"controlfreec_annotated_cn_autosomes.tsv.gz"))Parsed with column specification:
cols(
biospecimen_id = [31mcol_character()[39m,
status = [31mcol_character()[39m,
copy_number = [32mcol_double()[39m,
ploidy = [32mcol_double()[39m,
gene_symbol = [31mcol_character()[39m,
cytoband = [31mcol_character()[39m,
ensembl = [31mcol_character()[39m
)# read in the clinical file that has the disease type information
histologies_df <- readr::read_tsv(file.path("..", "..", "data",
"pbta-histologies.tsv"))Parsed with column specification:
cols(
.default = col_character(),
OS_days = [32mcol_double()[39m,
age_last_update_days = [32mcol_double()[39m,
normal_fraction = [32mcol_double()[39m,
tumor_fraction = [32mcol_double()[39m,
tumor_ploidy = [32mcol_double()[39m
)
See spec(...) for full column specifications.# get ATRT biospecimen identifiers
atrt_specimens <- histologies_df %>%
dplyr::filter(disease_type_new == "Atypical Teratoid Rhabdoid Tumor") %>%
dplyr::pull(Kids_First_Biospecimen_ID)# filter the gene-level CN information to ATRT samples and gene symbols
# containing 'SMARC'
cnvkit_atrt <- cnvkit_df %>%
dplyr::filter(biospecimen_id %in% atrt_specimens,
grepl("SMARCB1", gene_symbol))
cnvkit_atrttable(cnvkit_atrt$status)
loss
27 controlfreec_atrt <- controlfreec_df %>%
dplyr::filter(biospecimen_id %in% atrt_specimens,
grepl("SMARCB1", gene_symbol))
controlfreec_atrttable(controlfreec_atrt$status)
gain loss
1 22 Find the one sample that reports a gain in ControlFreeC and look at that sample in CNVkit.
controlfreec_gain_sample <- controlfreec_atrt %>%
dplyr::filter(status == "gain") %>%
dplyr::pull(biospecimen_id)
cnvkit_atrt %>% dplyr::filter(biospecimen_id == controlfreec_gain_sample)CNVkit marks this as a loss.
Check which samples are in both.
intersect(controlfreec_atrt$biospecimen_id, cnvkit_atrt$biospecimen_id) [1] "BS_0EN5R1AP" "BS_12ZB7R6A" "BS_1C4WCNQF" "BS_5WK2KXYW" "BS_6GN1FGTB" "BS_74A1TB03"
[7] "BS_8WJR4AA3" "BS_9KDA3HG6" "BS_A94A7EJ4" "BS_CKY2TTXG" "BS_DE26D072" "BS_GTAJFTTR"
[13] "BS_JDMM2XG4" "BS_NDSWS9TS" "BS_NXZSD2KK" "BS_P7CY04E7" "BS_QMHND3QK" "BS_R5G3QP11"
[19] "BS_SCQNN2A2" "BS_SWHQ8FDW" "BS_W23ERYW3" "BS_WCNPE0AC"length(unique(cnvkit_atrt$biospecimen_id))[1] 25length(unique(controlfreec_atrt$biospecimen_id))[1] 22