Skip to content

Commit

Permalink
Added correction for first questions of TP2
Browse files Browse the repository at this point in the history
  • Loading branch information
nalcala committed Nov 24, 2022
1 parent 12bcd67 commit a13fba8
Showing 1 changed file with 62 additions and 0 deletions.
62 changes: 62 additions & 0 deletions Practical2/Correction/Q1to4.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
library(tidyverse)
library(MOFA2)
# setup reticulate with correct python install
library(reticulate)
use_python(system("which python3.7",intern=T), required=TRUE)

# Q1
#a) read data
Data.Clin = read_tsv("medical_genomics_course/Practical2/Data/Data.Clin.txt")
RNA = read.table("medical_genomics_course/Practical2/Data/RNA.tsv",h=T,row.names=1)
DNAMeth_promoter = read.table("medical_genomics_course/Practical2/Data/DNAMeth_promoter.tsv",h=T,row.names=1)
DNAMeth_genebody = read.table("medical_genomics_course/Practical2/Data/DNAMeth_genebody.tsv",h=T,row.names=1)
DNAMeth_enhancer = read.table("medical_genomics_course/Practical2/Data/DNAMeth_enhancer.tsv",h=T,row.names=1)

#b)
head(Data.Clin)
head(RNA)
dim(RNA)
head(DNAMeth_promoter)
dim(DNAMeth_promoter)
head(DNAMeth_genebody)
dim(DNAMeth_genebody)
head(DNAMeth_enhancer)
dim(DNAMeth_enhancer)

#Q2
##a) densities
ggplot( data=tibble(value=unlist(RNA)), mapping = aes(x=value) ) + geom_histogram()
ggplot( data=tibble(value=unlist(DNAMeth_enhancer)), mapping = aes(x=value) ) + geom_histogram()
ggplot( data=tibble(value=unlist(DNAMeth_genebody)), mapping = aes(x=value) ) + geom_histogram()
ggplot( data=tibble(value=unlist(DNAMeth_promoter)), mapping = aes(x=value) ) + geom_histogram()
### not really Gaussian but not too widespread and globally symmetric
##b) create mofa object
mofa_untrained = create_mofa( list(RNA=as.matrix(RNA),Meth_pro=as.matrix(DNAMeth_promoter),Meth_bod=as.matrix(DNAMeth_genebody),Meth_enh=as.matrix(DNAMeth_enhancer)) )

#Q3
##a)
opt_train = get_default_training_options(mofa_untrained)
opt_train$convergence_mode ="slow"
opt_model = get_default_model_options(mofa_untrained)
opt_model$num_factors = 5

#b)
mofa_untrained = prepare_mofa(mofa_untrained,model_options = opt_model,training_options = opt_train)

#c)
mofa_trained = run_mofa(mofa_untrained)

#Q4
##a)
plot_variance_explained(mofa_trained)

##b)
plot_factors(mofa_trained)

##c)
plot_factor_cor(mofa_trained)

## clean session
rm(mofa_untrained)
rm(list = c("DNAMeth_enhancer", "DNAMeth_genebody", "DNAMeth_promoter","RNA","file","opt_model","opt_train"))
gc()

0 comments on commit a13fba8

Please sign in to comment.