library(lmQCM)
setwd(dirname(rstudioapi::getSourceEditorContext()$path))
options(stringsAsFactors = F)
dataset = {} # KIRP
########################################################################
#                      Clinical
########################################################################
dataset[['clinical']] = read.table('KIRP/clinical/nationwidechildrens.org_clinical_patient_kirp.txt', header = T, sep = '\t')
dataset[['clinical']] = dataset[['clinical']][3:dim(dataset[['clinical']])[1],]
colnames(dataset[['clinical']])
head(dataset[['clinical']])
View(dataset)
colnames(dataset[['clinical']])
dataset[['clinical']]$age_at_diagnosis
dataset[['clinical']]$age_at_diagnosis = strtoi(dataset[['clinical']]$age_at_diagnosis)
dataset[['clinical']]$last_contact_days_to = strtoi(dataset[['clinical']]$last_contact_days_to)
dataset[['clinical']]$death_days_to = strtoi(dataset[['clinical']]$death_days_to)
dataset[['clinical']]$tobacco_smoking_history_indicator = strtoi(dataset[['clinical']]$tobacco_smoking_history_indicator)
dataset[['clinical']]$tobacco_smoking_history_indicator
dataset[['clinical']]$age_at_diagnosis = strtoi(dataset[['clinical']]$age_at_diagnosis)
dataset[['clinical']]$last_contact_days_to = strtoi(dataset[['clinical']]$last_contact_days_to)
dataset[['clinical']]$death_days_to = strtoi(dataset[['clinical']]$death_days_to)
print('use valid \'death_days_to\' to replace \'last_contact_days_to\'')
# days_to_last_followup and days_to_death
dataset[['clinical']]$survival_days = dataset[['clinical']]$last_contact_days_to
dataset[['clinical']]$survival_days[!is.na(dataset[['clinical']]$death_days_to)] = dataset[['clinical']]$death_days_to[!is.na(dataset[['clinical']]$death_days_to)]
# extract useful columns
dataset[['clinical']] = dataset[['clinical']][, c('bcr_patient_barcode', 'gender','age_at_diagnosis', 'vital_status','survival_days')]
dataset[['clinical']] = dataset[['clinical']][complete.cases(dataset[['clinical']]),]
print(paste0('[clinical] ', dim(dataset[['clinical']])[1], ' complete rows found in clinical data.'))
dataset[['clinical']]
dataset[['TMB']] = read.table("KIRP/TMB/KIRP_TMB.csv", header = T, row.names = 1, sep = ",")
print(paste0('[TMB] ', dim(dataset[['TMB']])[1], ' complete rows found in clinical data.'))
########################################################################
#                      CNB
########################################################################
# Xiaohui Zhan:
# https://www.nature.com/articles/ng.3725
# a) For high quality CNVs ,the length of segmental region >20kb
#
# b) The number of probes spanning a CNV (a segmental region) to be
#    at least 10 to decrease false positives in calling CNVs.
#
# c) For a segmental region ,if the segment mean < |0.2|,this segmental
#    region should be discard.(Generally ,we using +/-0.2 as threshold
#    for a duplication/deletion. Because lots of noise will be introduced
#    from the microarray. The thresholds(+/- 0.2) were derived by examining
#    the distribution of segment mean values from tumor and normal samples)
SNP = read.table("LUAD/CNB/broad.mit.edu_PANCAN_Genome_Wide_SNP_6_whitelisted.seg", sep = "\t", header = T)
SNP$LENGTH = SNP$End - SNP$Start
SNP.filtered = SNP[(SNP$LENGTH >= 20000) &
(SNP$Num_Probes >= 10) &
(abs(SNP$Segment_Mean) >= 0.2),]
# SNP = read.table("data/UCSC Xena/broad.mit.edu_PANCAN_Genome_Wide_SNP_6_whitelisted.xena", sep = "\t", header = T)
SNP.filtered.sum = aggregate(SNP.filtered$LENGTH, by=list(Category=SNP.filtered$Sample), FUN=sum)
# 1Mb = 1,000 kb = 1,000,000 pb
SNP.filtered.sum$LENGTH_KB = SNP.filtered.sum$x/1000
# remove normal group 10A, 11A, ...
SNP.filtered.sum = SNP.filtered.sum[as.numeric(substr(SNP.filtered.sum$Category, 14, 15)) == 1, ] # only based on primary cancer (01A)
barcode = substr(SNP.filtered.sum$Category, 1, 12)
length(unique(barcode)) == length(barcode)
# samebarcode = names(table(barcode)[table(barcode)>=2])
# same = unlist(lapply(samebarcode, function(x) which(grepl(x, SNP.filtered.sum$Category))))
# SNP.filtered.sum$Category[same]
#### Get patients information
pinfo = read.table("LUAD/CNB/data/UCSC Xena/TCGA_phenotype_denseDataOnlyDownload.tsv", sep = "\t", header = T)
pinfo$barcode = substr(pinfo$sample, 1, 12)
SNP.filtered.sum$CANCER = pinfo$X_primary_disease[match(barcode, pinfo$barcode)]
study.abbr = read.table("KIRP/CNB/data/TCGA study abbreviations.tsv", sep = "\t", header = T)
study.abbr$Study.Name = tolower(study.abbr$Study.Name)
SNP.f2 = SNP.filtered.sum[SNP.filtered.sum$CANCER %in% study.abbr$Study.Name,]
SNP.f2$CANCER_ABBR = study.abbr$Study.Abbreviation[match(SNP.f2$CANCER, study.abbr$Study.Name)]
SNP.f3 = SNP.f2[SNP.f2$CANCER_ABBR %in% c("BLCA", "BRCA","CESC","HNSC","KIRC", "KIRP", "LIHC", "LUAD",
"LUSC", "OV", "PAAD","STAD"),]
table(SNP.f3$CANCER_ABBR)
dataset[['CNB']] = SNP.f3[SNP.f3$CANCER_ABBR == "KIRP",]
dataset[['CNB']] = data.frame(cbind(dataset[['CNB']]$Category, dataset[['CNB']]$LENGTH_KB))
colnames(dataset[['CNB']]) = c("barcode", "LENGTH_KB")
dataset[['CNB']]$barcode = unlist(lapply(dataset[['CNB']]$barcode, function(x) substr(x, 1, 12) ))
print(paste0('[CNB] ', dim(dataset[['CNB']])[1], ' complete rows found in clinical data.'))