#Copyright (c) 2016 Genome Research Ltd .
#Author : Tallulah Andrews <tallulandrews@gmail.com>
#This file is part of PoissonUMIs.
#PoissonUMIs is free software : you can redistribute it and/or modify it under
#the terms of the GNU General Public License as published by the Free Software
#Foundation; either version 2 of the License, or (at your option) any later
#version.
#This program is distributed in the hope that it will be useful, but WITHOUT
#ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#You should have received a copy of the GNU General Public License along with
#this program . If not , see <http://www.gnu.org/licenses/>.
PoisUMI_calc_discrepancy <- function(expr_mat) {
nrows = length(expr_mat[,1])
ncols = length(expr_mat[1,])
discrepancy = matrix(rep(0.0, times =nrows*ncols), nrow=nrows, ncol=ncols)
out <- .C("calc_discrepancy",expr_mat=as.integer(expr_mat), nc = as.integer(ncols), ng = as.integer(nrows), as.double(discrepancy));
output = matrix(out$discrepancy, nrow=nrows, ncol=ncols)
colnames(output) <- colnames(expr_mat)
rownames(output) <- rownames(expr_mat)
return(output)
}
PoisUMI_pearson_cor <- function(expr_mat, lambdas=NA) {
nrows = length(expr_mat[,1])
ncols = length(expr_mat[1,])
if (is.na(lambdas)) {
lambdas = PoisUMI_calc_lambdas(expr_mat);
}
cor_mat = matrix(rep(0.0, times = ncols*ncols), nrow=ncols, ncol=ncols)
out <- .C("calc_difference",expr_mat=as.integer(expr_mat), nc = as.integer(ncols), ng = as.integer(nrows), as.double(lambdas), as.double(cor_mat));
output = matrix(out$cor_mat, nrow=nrows, ncol=ncols)
colnames(output) <- colnames(expr_mat)
rownames(output) <- colnames(expr_mat)
return(output)
}
#data = read.table("/lustre/scratch108/compgen/team218/TA/Bergiers_Wafergen/Bergiers_Wafergen_Transcriptome_DirAdj_wMulti_FeatureCounts.txt", header=T)
#fit = PoisUMI_Fit_Full_Poisson(as.matrix(data))
#library(skellam)
#PoisUMI_distance(data[,1],data[,600],lambdas[,1],lambdas[,600])
#PoisUMI_distance(data_keep[2,],data_keep[1,],lambda_keep[2,],lambda_keep[1,])
# This is symetrical, D(1,2)=D(2,1), and in range [0,Inf], and D(1,1) == 0
# However this is slow for practical applications.
PoisUMI_distance <- function(Obs1, Obs2, Lambda1, Lambda2, alpha) {
Obs1 = alpha*Obs1
Obs2 = alpha*Obs2
# Calculate p-value for each pair of elements
difference = unlist(Obs1-Obs2);
thing = cbind(abs(difference), Lambda1, Lambda2)
density_of_tails <- function(val, l1, l2) {
if (val == 0) {
return(1)
} else {
p_upper = pskellam(floor(val),lambda1=l1,lambda2=l2, lower.tail=F)
p_lower = pskellam(-ceiling(val),lambda1=l1,lambda2=l2, lower.tail=T)
return(p_upper+p_lower)
}
}
ps = apply(thing,1, function(x){density_of_tails(x[1],l1=x[2],l2=x[3])})
# Combine p-values using Fisher's method
chi_sq = -2*sum(log(ps))
df = 2*length(ps)
pchisq(chi_sq, df=df, lower.tail=FALSE)
return(chi_sq/df)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.