Homec4science
Diffusion scATAC-seq (master) Browse

scATAC-seq/scripts/toy_data/analyse_data_sequence.R
411e7adaffa4master

Last Change
Raw File

analyse_data_sequence.R
View Options

setwd(file.path("/", "local", "groux", "scATAC-seq"))
# libraries
library(RColorBrewer)
library(motifStack)
library(TFBSTools)
# functions
source(file.path("scripts", "functions.R"))
#' Scores a sequence, in integer format, given
#' a probability matrix.
#' @pwm a probability matrix of dimensions :
#' 1st 4 for A, C, G, T
#' 2nd the number of positions
#' @param sequence a vector of integer encoding
#' the sequence A:1, C:2, G:3, T:4
#' @param prior a vector containing the prior
#' probabilities about each nucleotide. It may
#' be the base composition. It should sum up
#' to 1.
#' @return the sequence cumulative score.
#' @author Romain Groux
score.seq = function(sequence, pwm, prior=rep(0.25,4))
{
score = 0
l_seq = length(sequence)
l_pwm = ncol(pwm)
n_shifts = l_seq - l_pwm + 1
pwm = log(pwm)
prior = log(prior)
for(shift in 1:n_shifts)
{ score_ = 0
for(j in 1:l_pwm)
{ base = sequence[shift+j-1]
score_ = score_ + pwm[base, j] - prior[base]
}
score = score + exp(score_)
}
return(score)
}
#' Given a set of sequences and two motifs, this function
#' assigns to each sequence a class label.
#' The class label is determined by scoring a sequence
#' with both motifs. The best scoring motif determine
#' the class label.
#' @param sequences a matrix containing the sequences
#' in integer format. It has the following dimensions :
#' 1st the number of sequences
#' 2nd the sequence length
#' @param motif.1 a probability matrix containing the
#' 1st motif. It has the following dimensions :
#' 1st 4 for A, C, G, T
#' 2nd the number of positions
#' @param motif.2 a probability matrix containing the
#' 2st motif. It has the following dimensions :
#' 1st 4 for A, C, G, T
#' 2nd the number of positions
#' @return a list with 2 items
#' scores : a matrix containing the sequence scores. It
#' has the following dimensions :
#' 1st the number of sequences
#' 2nd 2 for each motif
#' labels : a vector containing the sequence labels.
#' @author Romain Groux
predict.class = function(sequences, motif.1, motif.2)
{ scores.1 = apply(sequences, 1, score.seq, motif.1)
scores.2 = apply(sequences, 1, score.seq, motif.2)
labels = rep(1, nrow(sequences))
labels[which(scores.1 < scores.2)] = 2
return(list("scores"=cbind(scores.1, scores.2), "labels"=labels))
}
#' Computes the x and y coordinates of the ROC curve
#' given the sequence scores and their corresponding
#' true class labels.
#' @param scores a matrix containing the scores for
#' each sequence (on the rows), for each class (on
#' the columns)
#' @param labels the true class labels for each
#' sequence.
#' @return a matrix with 2 columns containing the
#' x and y coordinates of the curve. The x coordinates
#' correspond to the specificity values and the y
#' coordinates the sensitivity values.
#' @author Romain Groux
get.roc.coord = function(scores, labels)
{ # number of points
n = length(labels)
# order by descending label for class 1 and then class2
ord = order(scores[,1], scores[,2], decreasing=T)
labels = labels[ord]
idx1 = which(labels == 1)
idx2 = which(labels == 2)
# true positive and negative discovery rates
true = vector(length=length(labels)) ; n_true = 0 ;
false = vector(length=length(labels)) ; n_false = 0 ;
for(i in 1:length(labels))
{ if(labels[i] == 1)
{ n_true = n_true + 1 }
else
{ n_false = n_false + 1 }
true[i] = n_true / (0.5*n) # 0.5 because as many class 1 as 2
false[i] = n_false / (0.5*n)
}
# to draw plot
m = cbind(true, false)
colnames(m) = c("y", "x")
return(m)
}
# seed
set.seed(20191007)
# path to the images for the logo
path.a = file.path("res/A.png")
path.c = file.path("res/C.png")
path.g = file.path("res/G.png")
path.t = file.path("res/T.png")
# 2000 sequences classified
n.seq = 2000
#number of time a classification was repeated
n.runs = 50
# number of classes searches
n.classes = 2
# the shifting freedom allowed
n.shifts = 90
# the expected dimensionality of the prob array to read
dim = c(n.seq, n.classes, n.shifts, 2) # 2000 seq, 2 classes, 90 shifts, 2 flips
# where the data are
dir.data = file.path("data",
"toy_data")
# where the results are
dir.results = file.path("results",
"toy_data")
# sequences
sequences = as.matrix(read.table(file.path(dir.data,
"simulated_sequences_2class_flip.mat"))) + 1
# true motifs
motif.true.1 = as.matrix(read.table(file.path(dir.data,
"simulated_sequences_2class_flip_motif1.txt")))
motif.true.2 = as.matrix(read.table(file.path(dir.data,
"simulated_sequences_2class_flip_motif2.txt")))
# true class labels
labels.true = as.matrix(read.table(file.path(dir.data,
"simulated_sequences_2class_flip_classes.txt")))[,1]
# AUC sequences with true motifs
scores.true = predict.class(sequences, motif.true.1, motif.true.2)
labels = scores.true$labels
auc.true = auc(labels.true, labels)
scores.true = scores.true$scores
# sequence scores with found motifs
scores.found = array(dim=c(n.runs, n.seq, 2), data=0)
# AUCs
auc.found = rep(0, n.runs)
# go over each run
for(i in 1:n.runs)
{ # get motifs found
file.motif = file.path(dir.results,
sprintf("simulated_sequences_2class_flip_class_model_%d.mat",
i))
motifs.found = read.sequence.models(file.motif)$models
# display logo
# par(mfrow=c(2,2))
# plot.logo(motifs.found[,,1], path.a, path.c, path.g, path.t, main=sprintf("run %d class1", i))
# plot.logo(motifs.found[,,2], path.a, path.c, path.g, path.t, main=sprintf("run %d class2", i))
# plot.logo(motif.true.1, path.a, path.c, path.g, path.t, main=sprintf("true class1", i))
# plot.logo(motif.true.2, path.a, path.c, path.g, path.t, main=sprintf("true class2", i))
# get AUC
# class 1 -> label 1, class 2 -> label2
scores1 = predict.class(sequences, motifs.found[,,1], motifs.found[,,2])
labels1 = scores1$labels
auc1 = auc(labels.true, labels1)
auc.found[i] = auc1
scores.found[i,,] = scores1$scores
# class 1 -> label 2, class 2 -> label1
# NO NEED ! AUC value is the same. This is simply an anti-classifier :)
# scores2 = predict.class(sequences, motifs.found[,,2], motifs.found[,,1])
# labels2 = scores2$labels
# auc2 = auc(labels.true, labels2)
# auc.found[i,2] = auc2
# scores.found[i,,2,] = scores2$scores
# coord1 = get.roc.coord(scores1$scores, labels.true)
# coord2 = get.roc.coord(scores2$scores, labels.true)
# coord.true = get.roc.coord(scores.true, labels.true)
# plot(coord1[,"x"], coord1[,"y"], lwd=2, type='l')
# lines(coord2[,"x"], coord2[,"y"], lwd=2, col="blue")
# lines(coord.true[,"x"], coord.true[,"y"], lwd=2, col="red")
# segments(0,0,1,1)
}
# # plot results
# X11(width=12, height=8)
#
# m = matrix(nrow=3, ncol=2,
# data=1:6, byrow=T)
#
# layout(m, heights=c(2,1,1))
#
# # boxplot AUCs
# boxplot(auc.found, ylim=c(0,1), main="AUCs", ylab="AUC",
# cex.lab=2, cex.main=2)
# abline(h=auc.true, col="red", lwd=2, lty=2)
#
# # plot all ROC
# coord.true = get.roc.coord(scores.true, labels.true)
# # empty plot
# plot(seq(0,1,by=0.2), seq(0,1,by=0.2), col=0,
# main="ROCs", xlab="specificity", ylab="sensitivity",
# cex.lab=2, cex.main=2)
# for(i in 1:n.runs)
# { # compute true positive and negative discovery rates
# coord.found = get.roc.coord(scores.found[i,,], labels.true)
# # plot found motifs ROC
# lines(coord.found[,"x"], coord.found[,"y"], lwd=0.5)
# }
# # plot diagonale line
# segments(0,0,1,1, lwd=2, lty=2)
# # plot true motifs ROC
# lines(coord.true[,"x"], coord.true[,"y"], lwd=3, col="red")
#
# # true motif class 1
# plot.logo(motif.true.1, path.a, path.c, path.g, path.t,
# main="True motif class 1", cex.main=2)
# # x-axis
# axis(1, at=1:ncol(motif.true.1), labels=1:ncol(motif.true.1))
# # x-axis
# axis(2, at=0:2, labels=0:2)
#
# # true motif class 2
# plot.logo(motif.true.2, path.a, path.c, path.g, path.t,
# main="True motif class 2", cex.main=2)
# # x-axis
# axis(1, at=1:ncol(motif.true.2), labels=1:ncol(motif.true.2))
# # x-axis
# axis(2, at=0:2, labels=0:2)
#
# # best motif found
# idx.best = which.max(auc.found)
# file.motif.best = file.path(dir.results,
# sprintf("simulated_sequences_2class_flip_class_model_%d.mat",
# idx.best))
# motif.found.best = read.sequence.models(file.motif.best)$models
# # best found motif class 1
# plot.logo(motif.found.best[,,1], path.a, path.c, path.g, path.t,
# main="Best found motif class 1", cex.main=2)
# # x-axis
# axis(1, at=1:ncol(motif.found.best[,,1]), labels=1:ncol(motif.found.best[,,1]))
# # x-axis
# axis(2, at=0:2, labels=0:2)
# # best found motif class 1
# plot.logo(motif.found.best[,,2], path.a, path.c, path.g, path.t,
# main="Best found motif class 2", cex.main=2)
# # x-axis
# axis(1, at=1:ncol(motif.found.best[,,2]), labels=1:ncol(motif.found.best[,,2]))
# # x-axis
# axis(2, at=0:2, labels=0:2)
#
# dev.off()
# plot results
png(filename=file.path(dir.results,
"simulated_sequences_2class_flip_auc_roc.png"),
units="in", res=720, width=12, height=6)
# X11(width=12, height=6)
par(mfrow=c(1,2),
mar=c(5.1,5.1,4.1,2.1))
# boxplot AUCs
boxplot(auc.found, ylim=c(0,1), main="AUC values", ylab="AUC",
cex.lab=2, cex.main=2)
abline(h=auc.true, col="red", lwd=2, lty=2)
# plot all ROC
coord.true = get.roc.coord(scores.true, labels.true)
# empty plot
plot(seq(0,1,by=0.2), seq(0,1,by=0.2), col=0,
main="ROC curves", xlab="specificity", ylab="sensitivity",
cex.lab=2, cex.main=2)
for(i in 1:n.runs)
{ # compute true positive and negative discovery rates
coord.found = get.roc.coord(scores.found[i,,], labels.true)
# plot found motifs ROC
lines(coord.found[,"x"], coord.found[,"y"], lwd=0.5)
}
# plot diagonale line
segments(0,0,1,1, lwd=2, lty=2)
# plot true motifs ROC
lines(coord.true[,"x"], coord.true[,"y"], lwd=3, col="red")
dev.off()
png(filename=file.path(dir.results,
"simulated_sequences_2class_flip_best_motifs.png"),
units="in", res=720, width=7, height=8)
# X11(width=7, height=10)
# best motif found
idx.best = which.max(auc.found)
file.motif.best = file.path(dir.results,
sprintf("simulated_sequences_2class_flip_class_model_%d.mat",
idx.best))
motif.found.best = read.sequence.models(file.motif.best)$models
pfm.found.best = get_pfm_list(motif.found.best, "Motif found")
rownames(motif.true.1) = rownames(motif.true.2) = c('A', 'C', 'G', 'T')
pfm.true.1 = new("pfm",
mat=motif.true.1,
name="True motif class 1")
pfm.true.2 = new("pfm",
mat=motif.true.2,
name="True motif class 2")
pfm.found.best.1 = new("pfm",
mat=motif.found.best[,,1],
name="Found motif class 1")
pfm.found.best.2 = new("pfm",
mat=motif.found.best[,,2],
name="Found motif class 2")
motifStack(c(pfm.true.1,
pfm.true.2,
pfm.found.best.1,
pfm.found.best.2),
layout="treeview")
dev.off()
c4science · Help