#***************************************************************************************************

#******************************* Make immunology FM **********************************************

#***************************************************************************************************

library(mclust)

library(data.table)

library(parallel)

source("/research/users/ppolonen/git_home/common_scripts/featurematrix/functions_generate_fm.R")

# WD

setwd("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/data/")

# GEXP

gexp=t(get(load("data9544_with_gene_symbols.RData")))

# annotation table

check=get(load("Hemap_immunology_Annotations_8304.Rdata"))

annot=read.delim("anno_coord_data9544_15pct_bw2.5_updated.txt", stringsAsFactors=F, header=T)

# listing files, not needed anymore

# f=list.files("/research/groups/sysgen/PROJECTS/HEMAP_IMMUNOLOGY/Annotations_immunology/", pattern = ".csv", full.names = T)

# f=f[-1]

# annot_normal=do.call(rbind, lapply(f, read.csv, header=F, stringsAsFactors=F, skip=1))

#******************************** filtering *************************************

exvivotreatments_allowed <- c("none", "na", "control", "activation", "differentiation", "differentiation followed by activation", " differentiation followed by LPS+IFNg", "differentiation followed by IFN", "differentiation followed by inflammatory cytokines", "differentiation followed by LPS", "differentiation followed by CD40L", "differentiation followed by Poly(I:C)", "differentiation (IL-4)", "differentiation with EPO", "anti-IgM", "IL-2", "IL-3", "IL-3+CpG")

rm_healthy = annot$Sample.type%in%c("NonCancerHealthy")&!(annot$In.vivo.treatment%in%c("na", "none", "no")&annot$Ex.vivo.treatment%in%exvivotreatments_allowed)

rm_cancer_prolif = annot$Sample.type%in%c("Cancer", "Prolif")&!annot$Ex.vivo.treatment=="none"

rm_celline = annot$Sample.type%in%c("CellLine")&!annot$Ex.vivo.treatment%in%c("none", "na", "control")

rm_treatments = annot$GSE.identifier..experiment.%in%c("GSE26661")|annot$GSM.identifier..sample.%in%c("GSM425497", "GSM425499")

rm_noncancer = !annot$Sample.type%in%c("Cancer", "Prolif", "CellLine", "NonCancerHealthy")

rm_mm_pbmc = annot$colorClass=="MM"&!annot$Sample.isolation=="CD138+ plasma cells"

# cluster exclusion

# AML,GSE7538,24

# CLL,GSE18866,GSE9250,32

# LP,GSE12453,GSE7345,6

# MM,GSE24147,GSE24522,18

# MP,GSE12079,GSE15811,13

# TCL,GSE14879,16 microdissected, do not remove

rm_cluster_differ=annot$GSE.identifier..experiment.%in%c("GSE7538", "GSE9250", "GSE12453", "GSE18866","GSE24522", "GSE7345", "GSE12079", "GSE15811", "GSE24147")&annot$Sample.type=="Cancer"

annot_left_out = annot

annot_left_out$reason_removed[rm_healthy]="Normal Cell sample, ex-vivo or in-vivo treated"

annot_left_out$reason_removed[rm_cluster_differ]="Cancer sample, outlier cluster"

annot_left_out$reason_removed[rm_celline]="Cell line sample, ex-vivo treated"

annot_left_out$reason_removed[rm_cancer_prolif|rm_treatments]="Cancer sample, ex-vivo treated"

annot_left_out$reason_removed[rm_noncancer]="NonCancer sample, not healthy"

annot_left_out$reason_removed[rm_mm_pbmc]="NonCancer sample, not healthy"

annot_left_out <- annot_left_out[(rm_noncancer|rm_healthy|rm_cancer_prolif|rm_celline|rm_treatments|rm_cluster_differ|rm_mm_pbmc),]

fwrite(annot_left_out, "hemap_1072_leftout_reasonremoved.txt", sep="\t")

annot <- annot[!(rm_healthy|rm_cancer_prolif|rm_celline|rm_treatments|rm_cluster_differ|rm_noncancer|rm_mm_pbmc),]

gexp <- gexp[,!(rm_healthy|rm_cancer_prolif|rm_celline|rm_treatments|rm_cluster_differ|rm_noncancer|rm_mm_pbmc)]

newSamples=annot[!annot[,1]%in%check[,1],]

fwrite(newSamples, "new_samples_in_hemap_8472.txt", sep="\t")

dim(check[check[,1]%in%annot[,1],])

#****************************************************************************************

# annot[rm_healthy&annot$GSM.identifier..sample.%in%annot_normal[,1],] # excluded from previous normals

# new

annot$CLASS[annot$Sample.type%in%c("NonCancerHealthy")]=gsub("NonCancerHealthy_|NonCancerHealthy|na_","", gsub("_StemCell_|Myeloid_|_na_na|1_na_na|2_na_na|Lymphoid_|_G.*.|_na_G0|_na_G1.*.|T-|B-|M1-Differentiating|Erythroid_","", annot$CLASS[annot$Sample.type%in%c("NonCancerHealthy")]))

annot$CLASS[annot$Sample.type%in%c("NonCancerHealthy")][annot$CLASS[annot$Sample.type%in%c("NonCancerHealthy")]=="na"]="LymphNode_GerminalCentre"

annot$CLASS[annot$CLASS=="CD4+Tcell"]="RestingCD4+Tcell"

annot$CLASS[annot$CLASS=="CD8+Tcell"]="RestingCD8+Tcell"

annot$MAINCLASS[annot$Sample.type%in%c("NonCancerHealthy")]="NonCancerHealthy"

annot$colorClass[annot$Sample.type%in%c("NonCancerHealthy")][annot$CLASS[annot$Sample.type%in%c("NonCancerHealthy")]=="na"]="Lymphoid"

# old

# annot$CLASS[match(annot_normal[,1], annot$GSM.identifier..sample.)]=annot_normal[,5]

# annot$MAINCLASS[match(annot_normal[,1], annot$GSM.identifier..sample.)]=annot_normal[,3]

# annot$colorClass[match(annot_normal[,1], annot$GSM.identifier..sample.)]=annot_normal[,4]

#*********************************** compute geometric mean ********************************

# get certain genes GEXP

rownames(gexp)=paste("N:GEXP:", rownames(gexp), sep="")

dat_a=gexp[grep("GZMA|PRF1|GNLY|GZMH|GZMM", rownames(gexp)),]

# dat_a=matrix[grep("GZMA|PRF1", rownames(matrix)),]

dat=2^dat_a+0.01

rownames(dat)

gm1=log2(t(apply(dat, 2, gm_mean)))

rownames(gm1)="CytolyticScore"

# also HLA

dat_a2=gexp[rownames(gexp)%in%c("N:GEXP:B2M", "N:GEXP:HLA-A", "N:GEXP:HLA-B", "N:GEXP:HLA-C"),]

dat2=2^dat_a2+0.01

rownames(dat2)

gm2=log2(t(apply(dat2, 2, gm_mean)))

rownames(gm2)="HLAIScore"

# also HLAII

dat_a3=gexp[rownames(gexp)%in%c("N:GEXP:HLA-DMA",

                "N:GEXP:HLA-DMB",

                "N:GEXP:HLA-DPA1",

                "N:GEXP:HLA-DPB1",

                "N:GEXP:HLA-DRA",

                "N:GEXP:HLA-DRB1"),]

dat3=2^dat_a3+0.01

rownames(dat3)

gm3=log2(t(apply(dat3, 2, gm_mean)))

rownames(gm3)="HLAIIScore"

classification1=data.frame(t(rep("medium", length(gm1))), stringsAsFactors = F)

zscore=as.numeric(scale(t(gm1)))

classification1[zscore>=1]="high"

classification1[zscore<=(-1)]="low"

rownames(classification1)="CytolyticScore" 

colnames(classification1)=colnames(gexp)

classification2=data.frame(t(rep("medium", length(gm2))), stringsAsFactors = F)

zscore=as.numeric(scale(t(gm2)))

classification2[zscore>=1]="high"

classification2[zscore<=(-1)]="low"

rownames(classification2)="HLAIScore" 

colnames(classification2)=colnames(gexp)

classification3=data.frame(t(rep("medium", length(gm3))), stringsAsFactors = F)

zscore=as.numeric(scale(t(gm3)))

classification3[zscore>=1]="high"

classification3[zscore<=(-1)]="low"

rownames(classification3)="HLAIIScore" 

colnames(classification3)=colnames(gexp)

classification=data.frame(t(rbind(classification1,classification2,classification3)), stringsAsFactors = F)

immunoscores=as.data.frame(t(rbind(gm1, gm2, gm3)))

immunoscoresfm=make.features(immunoscores, datatype="SAMP", prefix="")

colnames(immunoscoresfm)=colnames(gexp)

immunoscores_class_fm=make.features(classification, datatype="SAMP", prefix="")

colnames(immunoscores_class_fm)=colnames(gexp)

# excluding categorical here, they slow everything down!

l.data_list=list(gexp, immunoscoresfm, immunoscores_class_fm)

data_list=data.frame(do.call(rbind, l.data_list))

# ******************************** Infer cell fractions ***********************************

# cibersort

results=read.delim("CIBERSORT-Results.txt", row.names = 1, header=T, stringsAsFactors = F)

colnames(results)=paste0("N:SAMP:CIBERSORT_", gsub(" |-|\\.", "_", colnames(results)), "")

cibersort=t(results)

cibersort=cibersort[,colnames(cibersort)%in%colnames(gexp)]

MCP=get(load("MCP_counter_data.Rdata"))

rownames(MCP)=paste0("N:SAMP:", rownames(MCP), "")

MCP=MCP[,colnames(MCP)%in%colnames(cibersort)]

l.fractions=list(cibersort, MCP)

fractions=data.frame(do.call(rbind, l.fractions), stringsAsFactors = F)

#********************************** Clinical data ******************************

files=list.files(".", pattern=".info.tsv", full.names = T)

surv_data=do.call(rbind, lapply(files, read.delim, header=T, stringsAsFactors=F))

surv_data=surv_data[surv_data[,1]%in%colnames(gexp),]

surv_data=surv_data[!is.na(surv_data[,3]),]

surv_d=t(surv_data[,2:3])

colnames(surv_d)=surv_data[,1]

rownames(surv_d)=c("N:CLIN:OS_Time", "B:CLIN:OS_Status")

surv_d=surv_d[,match(colnames(gexp), colnames(surv_d))]

colnames(surv_d)=colnames(gexp)

#******************************* tumor percentage ********************************

Sys.setlocale(locale="C")

sorteds=read.delim("sorted_samples.txt", stringsAsFactors=F, header=F)

gsm=unlist(lapply(1:dim(sorteds)[1], function(i){

  annot$GSM.identifier..sample.[annot$Sample.isolation%in%sorteds[i,1]&annot$colorClass%in%sorteds[i,2]]

}))

add=annot$GSM.identifier..sample.[grepl("CD303", annot$Sample.isolation)]

sorted=t(annot$GSM.identifier..sample.%in%c(gsm, add))

rownames(sorted)="B:CLIN:CELLS_SORTED"

colnames(sorted)=colnames(gexp)

sorted[,grepl("Padiatr", annot$Additional.notes)]=1

# sorted[,annot$Additional.notes=="The leukemic blasts were sorted based on CD41, CD7, CD117, CD33, and CD34 antibodies as previously described (Klin. Padiatr. 217, 126-134)."] = 1

tumor_per=gsub("blast%: |>=|%|blast cell percentage: |t_cell_purity: |>|;|blast count, % of sample, -1=unavailable : ","", annot$Purity.Tumor.cell.content)

tumor_per[tumor_per=="high"]=80

tumor_per[as.numeric(tumor_per)>100]=100

tumor_per[tumor_per=="-1"]=0

tumor_per[tumor_per%in%c("n./a.", "na")]=NA

malt=read.delim("clinical_annotations_MALT.txt", stringsAsFactors=F, header=T)

replace=malt[match(colnames(gexp), malt$GSMID),]

tumor_per[grepl("Percentage of tumor", tumor_per)]=replace$X..Tumor[grepl("Percentage of tumor", tumor_per)]

tumor_percentage=data.matrix(t(as.numeric(tumor_per)))

rownames(tumor_percentage)="N:SAMP:BLAST_TUMOR_PERCENTAGE"

colnames(tumor_percentage)=colnames(gexp)

# T-cell percentages

T_per=data.matrix(t(as.numeric(replace$X..T.cells)))

rownames(T_per)="N:SAMP:TCELL_PERCENTAGE"

colnames(T_per)=colnames(gexp)

tissue_per=data.matrix(t(as.numeric(replace$X..Lung)))

rownames(tissue_per)="N:SAMP:TISSUE_PERCENTAGE"

colnames(tissue_per)=colnames(gexp)

#*************************** adding some lymphoma annotation *************************************

anno = read.delim("GSE10846_series_matrix_info_ipi_clean.txt", stringsAsFactors=F, header=T)

anno=anno[anno[,1]%in%annot[,1],]

anno=anno[match(annot[,1], anno[,1]),]

annot$In.vivo.treatment[annot[,1]%in%anno[,1]]=gsub("*.*: |;", "", anno$chemotherapy[annot[,1]%in%anno[,1]])

annot$In.vivo.treatment[annot$In.vivo.treatment%in%"NA"]=NA

annot$dlbcl_ipi=NA

annot$dlbcl_ipi[annot[,1]%in%anno[,1]]=anno$ipi[annot[,1]%in%anno[,1]]

CHOP=t(annot$In.vivo.treatment%in%"CHOP-Like Regimen"*1)

RCHOP=t(annot$In.vivo.treatment%in%"R-CHOP-Like Regimen"*1)

CHOP[is.na(annot$In.vivo.treatment)]=NA

RCHOP[is.na(annot$In.vivo.treatment)]=NA

rownames(CHOP)="B:CLIN:Chemotherapy_CHOP"

rownames(RCHOP)="B:CLIN:Chemotherapy_RCHOP"

colnames(CHOP)=colnames(gexp)

colnames(RCHOP)=colnames(gexp)

# add cytogenetic information

genetics_org=read.delim("AML_preBALL_cytogenetics_vectors.txt", stringsAsFactors=F, header=T, row.names=1)

genetics=genetics_org[rownames(genetics_org)%in%colnames(gexp),]

genetics=t(genetics)*1

rownames(genetics)=paste("B:CLIN:", "GENETICS_", rownames(genetics), "", sep="")

cytogenetic=annot$Cytogenetics

cytogenetic[cytogenetic%in%c("na", "n/a", "", " ")]=NA

cytogenetic[grepl("without|unknown|remainingcytogenetics|no del13q|crlf2 fish: Normal|crlf2 fish: n/a", cytogenetic)]=NA

cytogenetic[grepl("ormal", cytogenetic)]="normal_karyotype"

cytogenetic[grepl("MLL", cytogenetic)]="MLL"

cytogenetic=gsub("complex aberrant karyotype", "complex karyotype", cytogenetic)

cytogenetic=gsub("hyperdiploid karyotype", "hyperdiploid", cytogenetic)

cytogenetic=gsub("TAL$", "TAL1", cytogenetic)

cytogenetic=gsub("remaining cytogenetics|other abNormalities", "other", cytogenetic)

cytogenetic=gsub("fish:|trisomy 8 |;deletion|; complex karyotype|, complex karyotype|: positive| chromosomal aberrations|/API2-MALT1|/API2-MALT1 negative|deletion *.*: negative|/IGH-MALT1|, plus other| plus other", "", cytogenetic)

cytogenetic=gsub("trisomy ", "trisomy", cytogenetic)

cytogenetic=gsub("TEL deleted", "TEL_deleted", cytogenetic)

cytogenetic=gsub("p13.1", "p13", cytogenetic)

cytogenetic=gsub(";$", "", cytogenetic)

cytogenetic=gsub("\\+ ", "+", cytogenetic)

cytogenetic=gsub("i\\(", "inv(", cytogenetic)

cytogenetic=gsub("complex karyotype", "complex_karyotype", cytogenetic)

cytogenetic=gsub("^ ", "", cytogenetic)

annot$cytogenetic_clean=cytogenetic

cytogenetic_terms=sort(unique(unlist(strsplit(cytogenetic, " "))))

cytogenetic_terms=unlist(strsplit(cytogenetic, "; |; |  |, | |/"))

cytogenetic_terms=gsub(" ", "", cytogenetic_terms)

terms=table(cytogenetic_terms)

shared=names(terms)[terms>5]

cytogenetics=do.call(rbind, mclapply(shared,FIND_LOGICAL, cytogenetic, mc.cores=8))

colnames(cytogenetics)=colnames(gexp)

#***********************

# age annotations

age=gsub("*.*: |^ |;| yr| age| years|-.*.$| Years|d 32.8|d 54.8", "", annot$Age)

age[age%in%c("na", "n/a", "", "not available")]=NA

age[grepl("month|Month", age)]=signif(as.numeric(gsub(" months.*.| months| Months", "", age[grepl("month|Month", age)]))/12, 2)

age[age=="Adult"]=30

age[age=="Children"]=5

age[age=="pediatric"]=1

age=t(as.numeric(age))

age[age>100]=NA

rownames(age)="N:CLIN:AGE"

colnames(age)=colnames(gexp)

# gender annotations

gender=toupper(annot$Gender)

gender[gender%in%c("GENDER: NOT AVAILABLE", "GENDER: NA;", "SEX: UNKNOWN;")]=NA

gender=gsub(" |;", "", gender)

gender=gsub("GENDER:|SEX/AGE:|/.*.|SEX:", "", gender)

gender[gender%in%c("F", "FEMALE", "WOMAN")]="female"

gender[gender%in%c("M", "MALE", "MAN")]="male"

gender=t(gender)

rownames(gender)="C:CLIN:GENDER"

colnames(gender)=colnames(gexp)

# race annotations

race=toupper(annot$Race)

race[grep("AGE", race)]=NA

race[grep("OTHER", race)]="OTHER"

race[grep("AFRICAN|RACE: AA;|RACE: B;", race)]="AFRICAN"

race[grep("HISPANIC|RACE: H; ", race)]="HISPANIC"

race[grepl("EUROPEAN|CAUCASIAN|ANGLO-AMERICAN|WHITE|RACE: W;|RACE: C;", race)]="EUROPEAN"

race[grep("ASIAN", race)]="ASIAN"

race[!grepl("ASIAN|AFRICAN|HISPANIC|EUROPEAN|OTHER", race)]=NA

race=t(race)

rownames(race)="C:CLIN:RACE"

colnames(race)=colnames(gexp)

#*************************** adding some myeloma annotation *************************************

annomm = read.delim2("GSE24080_MM_clininfo_GSMid_clean.txt", stringsAsFactors=F, header=T)

annomm2 = read.delim2("GSE19784_MM_survival_GSMid_iss_clean.txt", stringsAsFactors=F, header=T)

annomm$b2m=gsub("<0.5", "0.5", annomm$b2m)

annomm$b2m=as.numeric(gsub(",", ".", annomm$b2m))

annomm$aspc=as.numeric(gsub(",", ".", annomm$aspc))

annomm$bmpc=as.numeric(gsub(",", ".", annomm$bmpc))

# first combine the two:

library(data.table)

combmm=data.frame(rbindlist(list(annomm, annomm2), fill = TRUE), stringsAsFactors = F)

combmm=combmm[match(colnames(gexp), combmm$accession),]

# now add these vectors to annot table

age[colnames(gexp)%in%combmm$accession]=signif(combmm$age[colnames(gexp)%in%combmm$accession], 3)

gender[colnames(gexp)%in%combmm$accession]=combmm$sex[colnames(gexp)%in%combmm$accession]

combmm$race[combmm$race%in%"other"]="OTHER"

combmm$race[combmm$race%in%"white"]="EUROPEAN"

race[colnames(gexp)%in%combmm$accession]=combmm$race[colnames(gexp)%in%combmm$accession]

# time and status:

surv_d[1,colnames(gexp)%in%combmm$accession]=signif(combmm$os_time[colnames(gexp)%in%combmm$accession], 3)

surv_d[2,colnames(gexp)%in%combmm$accession]=signif(combmm$os_censor[colnames(gexp)%in%combmm$accession], 3)

# pfs time and status

pfs=cbind(combmm$pfs_time, combmm$pfs_censor)

pfs=t(pfs)

colnames(pfs)=colnames(gexp)

rownames(pfs)=c("N:CLIN:PFS_Time", "B:CLIN:PFS_Status")

# other myeloma annotations:

otherMM_C=t(combmm[,c(4,8)])

rownames(otherMM_C)=paste0("C:CLIN:MM_", toupper(rownames(otherMM_C)), "")

colnames(otherMM_C)=colnames(gexp)

# numeric myeloma:

otherMM_N=t(combmm[,c(9:20)])

rownames(otherMM_N)=paste0("N:CLIN:MM_", toupper(rownames(otherMM_N)), "")

colnames(otherMM_N)=colnames(gexp)

add_mm=t(combmm$cytogenetic_abnormalities)

colnames(add_mm)=colnames(gexp)

rownames(add_mm)=c("B:CLIN:MM_CYTOGENETIC_ABNORMALITIES")

l.clin=list(sorted, surv_d, pfs, tumor_percentage,T_per,tissue_per, CHOP, RCHOP, gender, age,race, cytogenetics, genetics, add_mm, otherMM_N, otherMM_C)

clin=data.frame(do.call(rbind, l.clin), stringsAsFactors = F)

#****************************** make annotation clusters ******************************************

annot$acute=rep("other", dim(annot)[1])

annot$acute[annot$colorClass=="AML"|annot$colorClass=="pre-B-ALL"|annot$colorClass=="T-ALL"]="acute_leukemias"

annot$acute[annot$colorClass=="CLL"|annot$colorClass=="CML"]="chronic_leukemias"

annot$acute[grepl("NonCancer", annot$MAINCLASS)]="NonCancer"

annot$disease=rep("other", dim(annot)[1])

annot$disease[grepl("NonCancer", annot$Sub.maps.available)]="NonCancer"

annot$CLASS2=annot$CLASS

annot$CLASS2[grepl("CellLine_Myeloma", annot$CLASS2)]="CellLine_Myeloma"

annot$CLASS2[grepl("CellLine_Leukemia", annot$CLASS2)]="CellLine_Leukemia"

annot$CLASS2[grepl("CellLine_Lymphoma", annot$CLASS2)]="CellLine_Lymphoma"

annot$CLASS2[grepl("CellLine_mix", annot$CLASS2)]="CellLine_mix"

findthese=c("NonCancer", "Cancer_Leukemia", "Cancer_Myeloma", "Cancer_Lymphoma","CellLine_Leukemia","CellLine_Lymphoma","CellLine_Myeloma","Prolif_Lymphoproliferative_ALPS","Prolif_Lymphoproliferative_MPN", "Prolif_Myeloproliferative_LCH_LC", "Prolif_Myeloproliferative_MDS")

for(f in findthese){

  annot$disease[grepl(f, annot$MAINCLASS)]=f

}

annot$subclasses=rep("other", dim(annot)[1])

findthese=c("NonCancer", "Cancer_Myeloma", "Cancer_Lymphoma", "CellLine_Leukemia", "CellLine_Lymphoma","CellLine_Myeloma","Prolif_Lymphoproliferative_ALPS","Prolif_Lymphoproliferative_MPN", "Prolif_Myeloproliferative_LCH_LC", "Prolif_Myeloproliferative_MDS")

findthese2=c("T-ALL", "pre-B-ALL", "AML","CML","CLL", "BCL", "TCL", "B-Lymphoid", "T-Lymphoid","Lymphoid", "Myeloid", "Erythroid", "StemCell")

for(f in findthese){

  annot$subclasses[grepl(f, annot$MAINCLASS)]=f

}

for(f in findthese2){

  annot$subclasses[grepl(f, annot$colorClass)]=f

}

DLBCL=c("Cancer_Lymphoma_BCL_DLBCL_ABC", "Cancer_Lymphoma_BCL_DLBCL_GCB", "Cancer_Lymphoma_BCL_DLBCL_na")

annot$subclasses[annot$CLASS2%in%DLBCL]="BCL_DLBCL"

annot$CLASS=gsub("_na|_check|_testicular", "",annot$CLASS)

# lymphoma annotations

bLY=(1:nrow(annot)%in%grep("Lymphoma_BCL",annot$CLASS))&(!1:nrow(annot)%in%grep("CellLine",annot$CLASS))&(!1:nrow(annot)%in%grep("NonCancer",annot$CLASS))

tLY=(1:nrow(annot)%in%grep("Lymphoma_TCL",annot$CLASS))&(!1:nrow(annot)%in%grep("CellLine",annot$CLASS))&(!1:nrow(annot)%in%grep("NonCancer",annot$CLASS))

annot$CLASS=gsub("Cancer_", "", annot$CLASS)

annot$tbLY=NA

annot$tbLY[bLY|tLY]=annot$CLASS[bLY|tLY]

# these are the terms to look for

# table(annot$disease)

# table(annot$colorClass)

# table(annot$acute)

# table(annot$subclasses)

# table(annot$tbLY)

#*******************************************************************************************

#*************************** annotation clusters ********************************************

#****************************** make clusters ******************************************

# make immunological normal annotation vectors

annot$plotNormals = "Other"

lv=!annot$Sample.type%in%"NonCancerHealthy"

annot$plotNormals[lv]=""

lv=grepl("RestingBcell|NaiveBcell|MemoryBcell|BcellActivated", annot$CLASS)

annot$plotNormals[lv]="B cell"

lv=grepl("GerminalCentre", annot$CLASS)

annot$plotNormals[lv]="Germinal centre cell"

lv=grepl("PlasmaBcell", annot$CLASS)

annot$plotNormals[lv]="Plasma cell"

lv=grepl("Tcell|NaturalKillerCell", annot$CLASS)

annot$plotNormals[lv]="T/NK cell"

lv=grepl("DendriticCell", annot$CLASS)

annot$plotNormals[lv]="Dendritic cell"

lv=grepl("Langerhans", annot$CLASS)

annot$plotNormals[lv]="Langerhans cell"

lv=grepl("Eryth|Platelet", annot$CLASS)

annot$plotNormals[lv]="Erythroid"

lv=grepl("Monocyte", annot$CLASS)

annot$plotNormals[lv]="Monocyte"

lv=grepl("Macrophage", annot$CLASS)

annot$plotNormals[lv]="Macrophage"

lv=grepl("Neutrophil", annot$CLASS)

annot$plotNormals[lv]="Neutrophil"

lv=grepl("MyeloidProgenitor", annot$CLASS)

annot$plotNormals[lv]="Myeloid progenitor"

lv=grepl("HematopoieticStemCell", annot$CLASS)

annot$plotNormals[lv]="HSC"

lv=grepl("^Mononuclear", annot$CLASS)

annot$plotNormals[lv]="PBMC"

lv=grepl("LymphNode", annot$CLASS)

annot$plotNormals[lv]="Lymph node"

HLAplot_normals <- c("B cell", "Plasma cell", "T/NK cell", "Dendritic cell", "Erythroid", "Monocyte", "Macrophage", "Neutrophil", "Myeloid progenitor", "HSC")

cytolyticplot_normals <- c("PBMC", "Lymph node")

costimplot_normals <- c(HLAplot_normals, "PBMC", "Lymph node", "Langerhans cell", "Germinal centre cell")

annot$immunoNormals=annot$Category.specifying.lineage.tumor.origin

lv=grepl("CD8|CD8+TcellActivated", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="CD8+Tcell"

lv=grepl("NaturalKiller", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="NKCell"

lv=grepl("M2-Macrophage", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="M2-Macrophage"

lv=grepl("M1-Macrophage", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="M1-Macrophage"

lv=grepl("DendriticCell", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="DendriticCell"

lv=grepl("Monocyte", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="Monocyte"

lv=grepl("CD4+", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="CD4+Tcell"

lv=!annot$Sample.type%in%"NonCancerHealthy"

annot$immunoNormals[lv]=""

lv=grepl("Eryth|Platelet", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="Erythroid"

lv=grepl("CD3", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="Tcell"

lv=grepl("GerminalCentre", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="GerminalCentreCell"

lv=grepl("^Tcell$|^TcellActivated$|^TcellResting$", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="Tcell"

lv=grepl("^ActivatedBcell$|^RestingBcell$|^BcellActivated$", annot$Category.specifying.lineage.tumor.origin)

annot$immunoNormals[lv]="Bcell"

# annotated clusters

tbLY=FUN_MAKE_ALL(annot$tbLY, "annotated_class", annot$tbLY, 0)

subclasses=FUN_MAKE_ALL(annot$subclasses, "annotated_class", annot$subclasses, 0)

acute_chronic=FUN_MAKE_ALL(annot$acute, "annotated_class", annot$acute, 0)

colorClass=FUN_MAKE_ALL(annot$colorClass, "annotated_class", annot$colorClass, 0)

disease=FUN_MAKE_ALL(annot$disease, "annotated_class", annot$disease, 0)

tbLY=FUN_MAKE_ALL(annot$tbLY, "annotated_class", annot$tbLY, 0)

fullclass=FUN_MAKE_ALL(annot$CLASS2, "annotated_class", annot$CLASS2, 0)

immunoclass=FUN_MAKE_ALL(annot$immunoNormals, "annotated_class_immunoNormals", annot$immunoNormals, 0)

l.comparisons=list(disease, colorClass, acute_chronic,subclasses, tbLY, fullclass, immunoclass)

comparisons=do.call(rbind, l.comparisons)

comparisons=data.frame(data.matrix(comparisons[!duplicated(rownames(comparisons)),]), stringsAsFactors = F)

colnames(comparisons)=colnames(gexp)

# test if all rows are fine, should be >1 values

A=apply(comparisons, 1, unique)

B=unlist(lapply(A, function(d)sum(d%in%c(1,0))>=2))

if(!all(B))stop("Check comparisons, impossible comparisons made")

# categorical feats

class1=FUN_MAKE_CATEGORICAL(annot$tbLY, "annotated_class_BCL_TCL")

class2=FUN_MAKE_CATEGORICAL(annot$colorClass, "annotated_class_colorclass")

class3=FUN_MAKE_CATEGORICAL(annot$disease, "annotated_class_disease")

class4=FUN_MAKE_CATEGORICAL(annot$immunoNormals, "annotated_class_immunoNormals")

l.comparisons=list(class1, class2, class3, class4)

comparisons_cat=do.call(rbind, l.comparisons)

comparisons_cat=data.frame(comparisons_cat[!duplicated(rownames(comparisons_cat)),], stringsAsFactors = F)

colnames(comparisons_cat)=colnames(gexp)

#****************************************************************************************************

#******************** Next we start to create features of these individual data types ***************

#****************************************************************************************************

l.fm=list(data_list,clin,fractions, comparisons)

library(data.table)

fm=rbindlist(l.fm, use.names=F, fill=F)

fm=data.frame(fm, stringsAsFactors=F)

rownames(fm)=unlist(lapply(l.fm, rownames))

matrix=fm

# also add clinicaldata to annotations

numclin=t(data.matrix(clin[!grepl("^C:", rownames(clin)),]))

numchr=t(clin[grepl("^C:", rownames(clin)),])

colnames(numclin)=gsub(".:CLIN:|.:GEXP:|", "", colnames(numclin))

colnames(numchr)=gsub(".:CLIN:|.:GEXP:|", "", colnames(numchr))

annot_add=data.frame(numclin, numchr, stringsAsFactors = F)

fractions2=t(fractions)

colnames(fractions2)=gsub("N:SAMP:", "", colnames(fractions2))

annot2=data.frame(annot, annot_add, "CytolyticScore"=as.numeric(gm1) ,"HLAIScore"=as.numeric(gm2), "HLAIIScore"=as.numeric(gm3),classification,fractions2, stringsAsFactors = F)

# annot

clusters=read.delim("AML_15pct_BHSNE_mean-shift.txt", stringsAsFactors=F, header=T)

clusters=clusters[clusters$ID%in%annot$GSM.identifier..sample.,]

clusters_cancermap=clusters$X1.5..cluster

matrix_sub=matrix[,colnames(matrix)%in%clusters$ID]

annot_sub=annot[annot$GSM.identifier..sample.%in%clusters$ID,]

# TCGA clusters

cluster_mapping=read.delim("Table_TCGA_cluster_AML_cluster_assignment.txt", header=T, stringsAsFactors=F, sep="\t")

TCGA_cluster=rep("NA", dim(annot_sub)[1])

TCGA_cluster[clusters_cancermap%in%cluster_mapping[1,2]]="TCGA_AML_cluster_1"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[2,2]]="TCGA_AML_cluster_2"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[3:5,2]]="TCGA_AML_cluster_3"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[6:10,2]]="TCGA_AML_cluster_4"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[11,2]]="TCGA_AML_cluster_5"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[12:15,2]]="TCGA_AML_cluster_6"

TCGA_cluster[clusters_cancermap%in%cluster_mapping[16:17,2]]="TCGA_AML_cluster_7"

n=lapply(unique(TCGA_cluster), function(i)annot_sub$GSM.identifier..sample.[TCGA_cluster%in%i])

names(n)=unique(TCGA_cluster)

n=n[!unique(TCGA_cluster)%in%"NA"]

save(n, file="Hemap_immunology_TCGA_clusters.Rdata")

save(matrix, file="Hemap_immunology_fm.Rdata")

save(annot2, file="Hemap_immunology_Annotations.Rdata")

write.table(annot2,"Hemap_immunology_Annotations.tsv", sep="\t", col.names=T, row.names=F, quote=FALSE)

write.table(t(c("N:SAMP", as.character(colnames(matrix)))), file="Hemap_immunology_fm.tsv", sep="\t", col.names=F, row.names=F, quote=FALSE, append=F)

write.table(matrix, file="Hemap_immunology_fm.tsv", sep="\t", col.names=F, row.names=T, quote=FALSE, append=T)

# make a small fix here to harmonize survival data to months:

load("Hemap_immunology_Annotations.Rdata")

unique(cbind(annot2[!is.na(annot2$OS_Time), c(2,4)]))

annot2$OS_Time[!is.na(annot2$OS_Time)&annot2[,2]%in%c("GSE10846,GSE11318", "GSE10846", "GSE10846,GSE17372", "GSE11877")]=annot2$OS_Time[!is.na(annot2$OS_Time)&annot2[,2]%in%c("GSE10846,GSE11318", "GSE10846", "GSE10846,GSE17372", "GSE11877")]*12

# for myeloma, transform data to 5y survival to compare data sets:

modify=!is.na(annot2$OS_Time)&annot2[,2]%in%c("GSE16716,GSE24080")

find=annot2$OS_Time>60&modify

# change status to alive if dead later

find2=annot2$OS_Status==1&modify

annot2$OS_Status[find&find2]=0

annot2$OS_Time[find]=60

save(annot2, file="Hemap_immunology_Annotations.Rdata")

#****************************************************************************************

# This FM can then be used as a backbone for other FMs. GSVA and clusters must be added

#****************************************************************************************

matrix=get(load("Hemap_immunology_fm.Rdata"))

annot=get(load("Hemap_immunology_Annotations.Rdata"))

#********************************** Full map **********************************

clusters=read.delim("anno_coord_data9544_15pct_bw2.5_updated.txt", stringsAsFactors=F, header=T)

clusters=clusters[clusters$ID%in%annot$GSM.identifier..sample.,]

#*********************************** GSVA input ****************************

gsva=get(load("data8238_dufva_immunological_genes_updated_2016_GSVA_geneperm_lean_eFDR.Rdata"))

bindea=get(load("data8238_all_samples_dufva_bindea_2013_geneset_GSVA.Rdata"))

load("data8238_all_samples_Combined_pathway_signatures_210616_GSVA.Rdata")

gsva_es=rbind(gsva, bindea, gsva_es)

gsva_es=gsva_es[!duplicated(rownames(gsva_es)),]

# match cols gsva

gsva_es=data.frame(gsva_es[,match(colnames(matrix), colnames(gsva_es))])

colnames(gsva_es)=colnames(matrix)

#**************************** Cancermap clusters ****************************

clusters_cancermap=clusters$X2.5..cluster

cluster_cancermap=FUN_MAKE_ALL(clusters_cancermap, "cancermap_cluster", clusters_cancermap, 0)

subclasses=FUN_MAKE_ALL(annot$subclasses, "cancermap_cluster", clusters_cancermap, 0.8)

acute_chronic=FUN_MAKE_ALL(annot$acute, "cancermap_cluster", clusters_cancermap, 0.8)

colorClass=FUN_MAKE_ALL(annot$colorClass, "cancermap_cluster", clusters_cancermap, 0.8)

disease=FUN_MAKE_ALL(annot$disease, "cancermap_cluster", annot$disease, 0.8)

fullclass=FUN_MAKE_ALL(annot$CLASS2, "cancermap_cluster", annot$CLASS2, 0.8)

class_cancermap=FUN_MAKE_CATEGORICAL(clusters_cancermap, "cancermap_cluster")

l.comparisons=list(cluster_cancermap, subclasses, acute_chronic, colorClass, disease,fullclass, class_cancermap)

comparisons_cat=do.call(rbind, l.comparisons)

comparisons_cat=data.frame(comparisons_cat[!duplicated(rownames(comparisons_cat)),], stringsAsFactors = F)

colnames(comparisons_cat)=colnames(matrix)

# combine

l.fm=list(matrix, gsva_es, comparisons_cat)

library(data.table)

fm=rbindlist(l.fm, use.names=F, fill=F)

fm=data.frame(fm, stringsAsFactors=F)

rownames(fm)=unlist(lapply(l.fm, rownames))

# remove rows with few values or NAs

rm=apply(fm, 1, function(v)all(is.na(v)))

fm=fm[!rm,]

save(fm, file="Hemap_immunology_fm_cancermap.Rdata")

write.table(t(c("N:SAMP", as.character(colnames(fm)))), file="Hemap_immunology_fm_cancermap.tsv", sep="\t", col.names=F, row.names=F, quote=FALSE, append=F)

write.table(fm, file="Hemap_immunology_fm_cancermap.tsv", sep="\t", col.names=F, row.names=T, quote=FALSE, append=T)

#********************************** Lymphoma **********************************

matrix=get(load("Hemap_immunology_fm.Rdata"))

annot=get(load("Hemap_immunology_Annotations.Rdata"))

# annot

clusters=read.delim("Hemap_Lymphoma_15pct_genes_BHSNE_mean-shift.txt", stringsAsFactors=F, header=T)

clusters=clusters[clusters$ID%in%annot$GSM.identifier..sample.,]

matrix_sub=matrix[,colnames(matrix)%in%clusters$ID]

annot_sub=annot[annot$GSM.identifier..sample.%in%clusters$ID,]

load("data9544_LYMPHOMA_all_samples_Combined_pathway_drug_signatures_2017_GSVA.Rdata")

rownames(gsva_es)=gsub(" ", "_", rownames(gsva_es))

# match cols gsva

gsva_es=data.frame(gsva_es[,match(colnames(matrix_sub), colnames(gsva_es))])

colnames(gsva_es)=colnames(matrix_sub)

clusters_cancermap=clusters$X1.5..cluster

# comparisons

cluster_subtypes=FUN_MAKE_ALL(annot_sub$CLASS, "cancermap_cluster", clusters_cancermap, 0.8)

cluster_cancermap=FUN_MAKE_ALL(clusters_cancermap, "cancermap_cluster", clusters_cancermap, 0)

cluster_BCL_TCL=FUN_MAKE_ALL(annot_sub$tbLY, "cancermap_cluster", clusters_cancermap, 0.8)

class_cancermap=FUN_MAKE_CATEGORICAL(clusters_cancermap, "cancermap_cluster")

l.comparisons=list(cluster_subtypes, cluster_cancermap, cluster_BCL_TCL, class_cancermap)

comparisons_cat=do.call(rbind, l.comparisons)

comparisons_cat=data.frame(data.matrix(comparisons_cat[!duplicated(rownames(comparisons_cat)),]), stringsAsFactors = F)

# combine

l.fm=list(matrix_sub, data.frame(gsva_es), comparisons_cat)

library(data.table)

fm=rbindlist(l.fm, use.names=F, fill=F)

fm=data.frame(fm, stringsAsFactors=F)

rownames(fm)=unlist(lapply(l.fm, rownames))

# remove rows with few values or NAs

rm=apply(fm, 1, function(v)all(is.na(v)))

fm=fm[!rm,]

save(fm, file="Hemap_LYMPHOMA_immunology_fm.Rdata")

write.table(t(c("N:SAMP", as.character(colnames(fm)))), file="Hemap_LYMPHOMA_immunology_fm.tsv", sep="\t", col.names=F, row.names=F, quote=FALSE, append=F)

write.table(fm, file="Hemap_LYMPHOMA_immunology_fm.tsv", sep="\t", col.names=F, row.names=T, quote=FALSE, append=T)

#********************************** AML **********************************

matrix=get(load("Hemap_immunology_fm.Rdata"))

annot=get(load("Hemap_immunology_Annotations.Rdata"))

# annot

clusters=read.delim("AML_15pct_BHSNE_mean-shift.txt", stringsAsFactors=F, header=T)

clusters=clusters[clusters$ID%in%annot$GSM.identifier..sample.,]

clusters_cancermap=clusters$X1.5..cluster

matrix_sub=matrix[,colnames(matrix)%in%clusters$ID]

annot_sub=annot[annot$GSM.identifier..sample.%in%clusters$ID,]

# GSVA input

gsva=get(load("data9544_AML_all_samples_Combined_pathway_drug_signatures_2017_GSVA.Rdata"))

# match cols gsva

gsva_es=data.frame(gsva_es[,match(colnames(matrix_sub), colnames(gsva_es))])

colnames(gsva_es)=colnames(matrix_sub)

# TCGA clusters

cluster_mapping=read.delim("Table_TCGA_cluster_AML_cluster_assignment.txt", header=T, stringsAsFactors=F, sep="\t")

annot_sub$TCGA_cluster=rep("NA", dim(annot_sub)[1])

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[1,2]]="TCGA_AML_cluster_1"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[2,2]]="TCGA_AML_cluster_2"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[3:5,2]]="TCGA_AML_cluster_3"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[6:10,2]]="TCGA_AML_cluster_4"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[11,2]]="TCGA_AML_cluster_5"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[12:15,2]]="TCGA_AML_cluster_6"

annot_sub$TCGA_cluster[clusters_cancermap%in%cluster_mapping[16:17,2]]="TCGA_AML_cluster_7"

# comparisons

cluster_TCGA=FUN_MAKE_ALL(annot_sub$TCGA_cluster, "cancermap_cluster", annot_sub$TCGA_cluster, 0.9)

cluster_cancermap=FUN_MAKE_ALL(clusters_cancermap, "cancermap_cluster", clusters_cancermap, 0)

cluster_subtypes=FUN_MAKE_ALL(annot_sub$CLASS, "cancermap_cluster", clusters_cancermap, 0.9)

class_cancermap=FUN_MAKE_CATEGORICAL(clusters_cancermap, "cancermap_cluster")

class_TCGA=FUN_MAKE_CATEGORICAL(annot_sub$TCGA_cluster, "cancermap_cluster")

l.comparisons=list(cluster_TCGA, cluster_cancermap, cluster_subtypes, class_cancermap, class_TCGA)

comparisons_cat=do.call(rbind, l.comparisons)

comparisons_cat=data.frame(data.matrix(comparisons_cat[!duplicated(rownames(comparisons_cat)),]), stringsAsFactors = F)

# combine

l.fm=list(matrix_sub, data.frame(gsva_es), comparisons_cat)

library(data.table)

fm=rbindlist(l.fm, use.names=F, fill=F)

fm=data.frame(fm, stringsAsFactors=F)

rownames(fm)=unlist(lapply(l.fm, rownames))

# remove rows with few values or NAs

rm=apply(fm, 1, function(v)all(is.na(v)))

fm=fm[!rm,]

save(fm, file="Hemap_AML_immunology_fm.Rdata")

write.table(t(c("N:SAMP", as.character(colnames(fm)))), file="Hemap_AML_immunology_fm.tsv", sep="\t", col.names=F, row.names=F, quote=FALSE, append=F)

write.table(fm, file="Hemap_AML_immunology_fm.tsv", sep="\t", col.names=F, row.names=T, quote=FALSE, append=T)

# write.table(all_isolation, file="Hemap_all_isolation.tsv", sep="\t", col.names=F, row.names=F, quote=FALSE, append=F)