Mercurial > repos > vipints > deseq_hts

comparison deseq-hts_1.0/src/deseq-hts.sh @ 0:94a108763d9e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
deseq-hts version 1.0 wraps the DESeq 1.6.0
author vipints
date Wed, 09 May 2012 20:43:47 -0400
parents
children 8ab01cc29c4b
comparison
equal deleted inserted replaced
-1:000000000000 0:94a108763d9e
1 #/bin/bash
2 #
3 # This program is free software; you can redistribute it and/or modify
4 # it under the terms of the GNU General Public License as published by
5 # the Free Software Foundation; either version 3 of the License, or
6 # (at your option) any later version.
7 #
8 # Copyright (C) 2009-2012 Max Planck Society
9 #
10
11 set -e
12
13 PROG=`basename $0`
14 DIR=`dirname $0`
15
16 . ${DIR}/../bin/deseq_config.sh
17
18 echo
19 echo ${PROG}: FML http://galaxy.fml.mpg.de Galaxy wrapper for the DESeq version $DESEQ_VERSION.
20 echo
21 echo DESeq performs differential expression testing from RNA-Seq measurements.
22 echo
23
24 ANNO_INPUT=${1}
25 shift
26 DESEQ_RES_FILE=${1}
27 shift
28 GENES_FN=${1}
29 shift
30
31 mkdir -p `dirname $GENES_FN`
32
33 echo %%%%%%%%%%%%%%%%%%%%%%%
34 echo % 1. Data preparation %
35 echo %%%%%%%%%%%%%%%%%%%%%%%
36 echo
37 echo load the genome annotation in GFF3 format and create an annotation object
38 echo
39 export PYTHONPATH=$PYTHONPATH:${SCIPY_PATH}
40 ${PYTHON_PATH} ${DIR}/../tools/ParseGFF.py ${ANNO_INPUT} ${GENES_FN}
41 ${DIR}/../bin/genes_cell2struct ${GENES_FN} 2>&1
42 echo
43 echo genome annotation stored in $GENES_FN
44
45 echo
46 echo %%%%%%%%%%%%%%%%%%%%
47 echo % 2. Read counting %
48 echo %%%%%%%%%%%%%%%%%%%%
49 echo
50
51 echo counting reads overlapping exons using given alignments
52 for REPLICATE_GROUP in $@
53 do
54 IFS=':'
55 for BAM_FILE in ${REPLICATE_GROUP}
56 do
57 echo
58 if [ ! -f ${BAM_FILE}.bai ]
59 then
60 echo "Indexing $BAM_FILE"
61 ${SAMTOOLS_DIR}/samtools index $BAM_FILE
62 else
63 echo "$BAM_FILE already indexed"
64 fi
65 echo
66 done
67 done
68 tmpfile=`mktemp --tmpdir=/tmp`
69
70 echo "${DIR}/../bin/get_read_counts ${GENES_FN} $tmpfile $@"
71 ${DIR}/../bin/get_read_counts ${GENES_FN} $tmpfile $@ 2>&1
72
73 echo
74 echo %%%%%%%%%%%%%%%%%%%%%%%%%%%
75 echo % 3. Differential testing %
76 echo %%%%%%%%%%%%%%%%%%%%%%%%%%%
77 echo
78
79 echo testing genes for differential expression using given alignments
80
81 echo "cat ${DIR}/../src/difftest_deseq.R | $R_PATH --slave --args $tmpfile ${DESEQ_RES_FILE} $#"
82 cat ${DIR}/../src/difftest_deseq.R | $R_PATH --slave --args $tmpfile ${DESEQ_RES_FILE} $# 2> /dev/null
83
84 rm $tmpfile ${tmpfile}_COUNTS.tab ${tmpfile}_CONDITIONS.tab
85 echo
86 echo %%%%%%%%
87 echo % Done %
88 echo %%%%%%%%
89 echo