preprocessing.R

#190814 colData 
#remove columns with 1 factor 
#preprocessing
#check if multiple colnames after shortening are there 
countBRCA <- readRDS("GDCdata/countBRCA.rds")
methylation <- readRDS("P:/mmRmeta/GDCdata/processed_data_brca/tcga-methyl-brca.rds")
methylation <- filter_duplicates(methylation) %>% filter_zero
#preprocess your the assay in SummarizedExperiment, in this case gene expression is in assay one. This function concists of 2 parts:
#filter_duplicated: filters out duplicated gene names
#filter_rare 
subBRCA <- preprocess_assay(countBRCA)
#merge expression of duplicated patients/samples
subBRCA <- merge_duplicated(subBRCA, assay_index = 1, method = "median", shorten_rownames = T, original_names = T, parallel = T )
#preprocess coldata / metadata of SE, consits of:
#factorize_columns: make factors out of columns 
#filter out columns with only NA val
#filter out columns with only 1 factor level
subBRCA <- preprocess_coldata(subBRCA)
#keep only samples that are Primary Solid Tumors and female 
subBRCA <- filter_for_factor(subBRCA, column_value = "Primary solid Tumor") %>% filter_for_factor(. , column_value = "female")
#reorder factor levels by counts
subBRCA <- reorder_all(subBRCA)
#preprocess subset again
subBRCA <- preprocess_assay(subBRCA)
#use siber on SE
siberBRCA <- use_siber(subBRCA, 1) 
#rearange out of siber
siberBRCA <- siberBRCA %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()
#get names of genes with siber >= 1.5
geneBRCA <- rownames(siberBRCA[siberBRCA$BI >= 1.5,])
#subset SE by names, add normalized assay, log transform assay
subBRCA <- subBRCA[geneBRCA,] %>% add_assay_norm(., 1, add_info_metadata = F) %>% transform_assay(.)
#use mclust and verify output
clust <- use_mclust(subBRCA, 2, G = c(1:2), modelNames = "V")
clust <- verify_mclust(clust)
clust <- verify_mclust(clust, min_proportion = 0.05)
#subset data by names in clust
subBRCA <- filter_fold_change(clust, subBRCA)
#add classification
subBRCA <- add_assay_classification(subBRCA, clust)


#21.08.19, for every PAM subtype
lumACarcinoma <- filter_for_factor(subBRCA, column_value = "LumA" ) %>% preprocess_assay()
LumBCarcinoma <- filter_for_factor(subBRCA, column_value = "LumB") %>% preprocess_assay()
BasalCarcinoma <- filter_for_factor(subBRCA, column_value = "Basal") %>% preprocess_assay()
Her2Carcinoma <- filter_for_factor(subBRCA, column_value = "Her2") %>% preprocess_assay()
NormalCarcinoma <- filter_for_factor(subBRCA, column_value = "Normal") %>% preprocess_assay()

siberlumA <- use_siber(lumACarcinoma, 1)
siberlumA <- siberlumA %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()
siberLumB<- use_siber(LumBCarcinoma, 1)
siberLumB <- siberLumB %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()
siberBasal <- use_siber(BasalCarcinoma, 1)
siberBasal <- siberBasal %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()
siberHer2<- use_siber(Her2Carcinoma, 1)
siberHer2 <- siberHer2 %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()
siberNormal<- use_siber(NormalCarcinoma, 1)
siberNormal<- siberNormal %>% tibble::rownames_to_column() %>% arrange(., -BI,) %>% tibble::column_to_rownames()

genelumA <- rownames(siberlumA[siberlumA$BI >= 1.5,])
genelumB <- rownames(siberLumB[siberLumB$BI >= 1.5,])
geneBasal <- rownames(siberBasal[siberBasal$BI >= 1.5,])
geneHer2 <- rownames(siberHer2[siberHer2$BI >= 1.5,])
geneNormal <- rownames(siberNormal[siberNormal$BI >= 1.5,])

#extract genes over 1.5
lumACarcinoma <- lumACarcinoma[genelumA,] %>% add_assay_norm(., 1, add_info_metadata = F) %>% transform_assay(.)
LumBCarcinoma <- LumBCarcinoma[genelumB,] %>% add_assay_norm(., 1, add_info_metadata = F)  %>% transform_assay(.)
BasalCarcinoma <- BasalCarcinoma[geneBasal,] %>% add_assay_norm(., 1, add_info_metadata = F)  %>% transform_assay(.)
Her2Carcinoma <- Her2Carcinoma[geneHer2,] %>% add_assay_norm(., 1, add_info_metadata = F)  %>% transform_assay(.)
NormalCarcinoma <- NormalCarcinoma[geneNormal,] %>% add_assay_norm(., 1, add_info_metadata = F)  %>% transform_assay(.)

lumAClust <- use_mclust(lumACarcinoma, 2, parallel = T, G = c(1:2), modelNames = "V")
lumAClust <- verify_mclust(lumAClust)
lumAClust <- verify_mclust(lumAClust, min_proportion = 0.05)
lumAClust <- verify_mclust(lumAClust)
lumAClust <- filter_unimodal(lumAClust)
subA <- filter_fold_change(lumAClust, summarized_experiment = lumACarcinoma)

exp <-  assay(subBRCA, 2)["ENSG00000256618",]
meta <- unnest_dataframe(colData(subBRCA))
meta$exp <- exp


########## GO Plot install.packages('GOplot')
install.packages('GOplot')
plot_heatmap(subBRCA, 2, col_annotaion = a)

coldat <- unnest_dataframe(colData(subBRCA))
anno <- calculate_annotation(unnest_dataframe(colData(subBRCA)), columns = "subtype_BRCA_Subtype_PAM50")
expression <- assay(subBRCA, 2)
a <- ComplexHeatmap::HeatmapAnnotation(select(coldat, c("subtype_BRCA_Subtype_PAM50", "primary_diagnosis")))
ComplexHeatmap::Heatmap(matrix = expression, top_annotation = a, show_column_names = F, show_row_names = F, show_column_dend = F)
a <- methylation[1:500,]
a <- assay(a, 1)
a <- t(a)
b <- data.frame(a) %>% dplyr::select_if(~!all(is.na(.)))
plot_heatmap(a, 1, col_annotation = "subtype_BRCA_Subtype_PAM50", show_column_names = F, show_row_names = F, show_row_dend = F, show_column_dend = F)

a <- t(assay(subBRCA, 2)[1:100,1:100])
b <- dist(scale(a), method = "euclidean")
c <- hclust(b, method = "ward.D2")

a <- add_coldata_classification(subBRCA, 3, T)
b <- unnest_dataframe(colData(a)) %>% group_by()