Commit 6c932ed0 authored by khalid's avatar khalid
Browse files

Disbled B2, Salmon and Hapmap indexing

parent 90e52ae7
......@@ -58,11 +58,11 @@ if [ $RECAT == "no" ]; then
### index contigs for the selected tool
case "$TOOL" in
B) toolidx="ALL_transcripts_bowtie_index"; if [ ! -d $out/$toolidx ]; then mkdir $out/$toolidx; bowtie-build --offrate 3 $out/ALL_transcripts.fasta $out/$toolidx/$toolidx ;fi ;;
B2) toolidx="ALL_transcripts_bowtie2_index"; if [ ! -d $out/$toolidx ]; then mkdir $out/$toolidx; bowtie2-build --offrate 3 $out/ALL_transcripts.fasta $out/$toolidx/$toolidx ;fi ;;
#B2) toolidx="ALL_transcripts_bowtie2_index"; if [ ! -d $out/$toolidx ]; then mkdir $out/$toolidx; bowtie2-build --offrate 3 $out/ALL_transcripts.fasta $out/$toolidx/$toolidx ;fi ;;
K) toolidx="ALL_transcripts_kallisto_index"; if [ ! -f $out/$toolidx ]; then kallisto index -i $out/$toolidx $out/ALL_transcripts.fasta; fi ;;
S) toolidx="ALL_transcripts_salmon_index" ; if [ ! -d $out/$toolidx ]; then salmon --no-version-check index -t $out/ALL_transcripts.fasta -i $out/$toolidx; fi ;;
#S) toolidx="ALL_transcripts_salmon_index" ; if [ ! -d $out/$toolidx ]; then salmon --no-version-check index -t $out/ALL_transcripts.fasta -i $out/$toolidx; fi ;;
R) toolidx="ALL_transcripts_rapmap_index" ; if [ ! -d $out/$toolidx ]; then rapmap quasiindex -t $out/ALL_transcripts.fasta -p -x $PROCESSORS -i $out/$toolidx; fi ;;
H) toolidx="ALL_transcripts_hpg_index" ; if [ ! -d $out/$toolidx ]; then mkdir $out/$toolidx; hpg-aligner build-sa-index -g $out/ALL_transcripts.fasta -i $out/$toolidx; fi ;;
#H) toolidx="ALL_transcripts_hpg_index" ; if [ ! -d $out/$toolidx ]; then mkdir $out/$toolidx; hpg-aligner build-sa-index -g $out/ALL_transcripts.fasta -i $out/$toolidx; fi ;;
esac
echo -e "\nIndex built : $out/$toolidx\n"
......
#Rscript --vanilla calc_edgeR_TMM.R file.txt
calcFactorWeighted <- function (obs, ref, libsize.obs = NULL, libsize.ref = NULL, logratioTrim = 0.3,
sumTrim = 0.05, doWeighting = TRUE, Acutoff = -1e+10)
{
if (all(obs == ref))
return(1)
obs <- as.numeric(obs)
ref <- as.numeric(ref)
if (is.null(libsize.obs))
nO <- sum(obs)
else nO <- libsize.obs
if (is.null(libsize.ref))
nR <- sum(ref)
else nR <- libsize.ref
logR <- log2((obs/nO)/(ref/nR))
absE <- (log2(obs/nO) + log2(ref/nR))/2
v <- (nO - obs)/nO/obs + (nR - ref)/nR/ref
fin <- is.finite(logR) & is.finite(absE) & (absE > Acutoff)
logR <- logR[fin]
absE <- absE[fin]
v <- v[fin]
n <- sum(fin)
loL <- floor(n * logratioTrim) + 1
hiL <- n + 1 - loL
loS <- floor(n * sumTrim) + 1
hiS <- n + 1 - loS
keep <- (rank(logR) >= loL & rank(logR) <= hiL) & (rank(absE) >=
loS & rank(absE) <= hiS)
if (doWeighting)
2^(sum(logR[keep]/v[keep], na.rm = TRUE)/sum(1/v[keep],
na.rm = TRUE))
else 2^(mean(logR[keep], na.rm = TRUE))
}
calcFactorQuantile <- function (data, lib.size, p = 0.75)
{
y <- t(t(data)/lib.size)
f <- apply(y, 2, function(x) quantile(x, p = p))
}
calcNormFactors <- function (object, method = c("TMM", "RLE", "upperquartile", "none"),
refColumn = NULL, logratioTrim = 0.3, sumTrim = 0.05, doWeighting = TRUE,
Acutoff = -1e+10, p = 0.75)
{
if (is(object, "DGEList")) {
x <- as.matrix(object$counts)
lib.size <- object$samples$lib.size
}
else {
x <- as.matrix(object)
lib.size <- colSums(x)
}
method <- match.arg(method)
allzero <- rowSums(x > 0) == 0
if (any(allzero))
x <- x[!allzero, , drop = FALSE]
if (nrow(x) == 0 || ncol(x) == 1)
method = "none"
f <- switch(method, TMM = {
f75 <- calcFactorQuantile(data = x, lib.size = lib.size,
p = 0.75)
if (is.null(refColumn)) refColumn <- which.min(abs(f75 -
mean(f75)))
if (length(refColumn) == 0 | refColumn < 1 | refColumn >
ncol(x)) refColumn <- 1
f <- rep(NA, ncol(x))
for (i in 1:ncol(x)) f[i] <- calcFactorWeighted(obs = x[,
i], ref = x[, refColumn], libsize.obs = lib.size[i],
libsize.ref = lib.size[refColumn], logratioTrim = logratioTrim,
sumTrim = sumTrim, doWeighting = doWeighting, Acutoff = Acutoff)
f
}, RLE = calcFactorRLE(x)/lib.size, upperquartile = calcFactorQuantile(x,
lib.size, p = p), none = rep(1, ncol(x)))
f <- f/exp(mean(log(f)))
if (is(object, "DGEList")) {
object$samples$norm.factors <- f
return(object)
}
else {
return(f)
}
}
args = commandArgs(trailingOnly=TRUE)
z = read.table(args[2],sep=",",header=T)
calcNormFactors(z)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment