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.'))