--- a
+++ b/data/.Rhistory
@@ -0,0 +1,80 @@
+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.'))