n_integration

Multi-omics data integration for multiple measurements

from the same individuals. In data matrices, samples are

in columns and features in rows.

Datasets used in the study from SSC available in GEO

under accessions: GSE23682 (copynumber),

GSE37772 (expression) and GSE27044 (methylation)

Samplenames with xxxxxx.p or xxxxxx.s to indicate

autist proband and their unaffected sibling.

Datasets preprocessed and normalized accordngly

for each datatype prior to this point.

Install packages

x=c("mixOmics","org.Hs.eg.db","GenomicRanges","annotatr","AnnotationDbi","ReactomePA","clusterProfiler","enrichplot")
lapply(x,BiocManager::install(x),character.only=T)
lapply(x,library,character.only=T)

Prepare the data, for expression and methylation select 2000 most

variable genes / probes. One option is to run differential analysis

for each datatype and select N most differentially expressed genes / probes.

var=rowVars(exp)
ord=order(var, decreasing=T)
exp=exp[ord[1:2000],]
var=rowVars(meth)
ord=order(var, decreasing=T)
meth=meth[ord[1:2000],]

X = list("Expression"=exp, "Methtylation"=meth, "Copynumber"=cn)
Y = as.factor(substr(colnames(exp), 7,7))
levels(Y)=c("ctrl","asd")
names(Y)=colnames(exp)

Prepare pairs of datatypes for initial analysis to detect correlation

levels between datatypes

x1=list(t(X$Expression),t(X$Copynumber))
names(x1)=c("Expression", "Copynumber")
list.keepx1 <- list(Expression = c(20,20), Copynumber=c(5,5))
x2=list(t(X$Expression),t(X$Methylation))
names(x2)=c("Expression","Methylation")
list.keepx2 <- list(Expression = c(20,20), Methylation=c(20,20))
x2=list(t(X$Methylation),t(X$Copynumber))
names(x3)=c("Methylation","Copynumber")
list.keepx3 <- list(Methylation = c(20,20), Copynumber=c(5,5))

Perform block.plsda for each pair of data to see the level of correlation

blockpls1=block.plsda(x1, Y, ncomp=2)
plotDiablo(blockpls1)
perf1=perf(blockpls1, validation = 'Mfold', folds=3, nrepeat = 10)
plot(perf1,legend.position="horizontal")

blockpls2=block.plsda(x2, Y, ncomp=2)
plotDiablo(blockpls2)
perf2=perf(blockpls2, validation = 'Mfold', folds=3, nrepeat = 10)
plot(perf2, legend.position="horizontal", col=color.mixo(5:7))

blockpls3=block.plsda(x3, Y, ncomp=2)
plotDiablo(blockpls3)
perf3=perf(blockpls3, validation = 'Mfold', folds=3, nrepeat = 10)
plot(perf3, legend.position="horizontal", col=color.mixo(5:7))

Select optimal number of features and components (generally components = N(response)-1)

keep_exp=tune(method="splsda", X=t(exp),Y=Y, ncomp=2, test.keepX=c(5,50,100), nrepeat=5)
plot(keep_exp)
keep_m=tune(method="splsda", X=t(meth),Y=Y, ncomp=2, test.keepX=c(5,50,100), nrepeat=5)
plot(keep_m)
keep_cn=tune(method="splsda", X=t(cn),Y=Y, ncomp=2, test.keepX=c(5,50,100), nrepeat=5)
plot(keep_cn)

List the optimal number of features to keep from each component

list.keepx <- list(Expression=keep_exp$choice.keepX,Methylation = keep_m$choice.keepX, keep_cn$choice.keepX)

Design matrices for block.splsda, select correlation based on the block.plsda results

design_x1=matrix(c(0,0.3,0.5,0.3,0,0.5,0.5,0.5,0), nrow=3, ncol=3, dimnames=list(c(names(x1),"Y"), c(names(x1),"Y")))

splsda1=block.splsda(X, Y, keepX=list.keepx, ncomp=2, design = design_x1)

Plot results

plotDiablo(splsda1)
circosPlot(splsda1, cutoff=0.5, size.variables = 0.4)
cimDiablo(splsda1)

Selected variables

variables_x1=selectVar(splsda1, comp=1)
variables_x2=selectVar(splsda1, comp=2)

Make a granges object of the selected copynumber segments

sel_cn=c(variables_x1$Copynumber$name,variables_x2$Copynumber$name)
chrom=unlist(lapply(strsplit(sel_cn, ":"),[,1))
pos=unlist(lapply(strsplit(sel_cn, ":"),[,2))
start=unlist(lapply(strsplit(pos, "-"),[,1))
end=unlist(lapply(strsplit(pos, "-"),[,2))
selected_gr=GRanges(seqnames=chrom, ranges=IRanges(start=as.numeric(start), end=as.numeric(end)), strand=rep("*",length(chrom)))

Annotate segments

annots=c("hg38_cpgs", "hg38_basicgenes")
annotations = build_annotations(genome = 'hg38', annotations = annots)
cn_ann=annotate_regions(regions=selected_gr, annotations=annotations, ignore.strand=T, quiet=FALSE, minoverlap=1)
summarize_annotations(cn_ann)
annot_df=as.data.frame(cn_ann$annot)
genes_cn=unique(annot_df$symbol)
genes_cn=genes_cn[!is.na(genes_cn)]

List of all selected features as genes

genelist=c(variables_x1$Expression$name, variables_x2$Expression$name,
variables_x1$Methylation$name, variables_x2$Methylation$name,
genes_cn)

Pathway analysis using ReactomePA

data(egSymb)
entrez_genes=sym2eg(genelist)
res=enrichPathway(gene=entrez_ids, organism="human", pvalueCutoff=0.1, qvalueCutoff=0.1,readable=T)

Plot results

dotplot(res)