added technical replicate merging to dge_utils

dge_workflow/README.md

@@ -23,8 +23,8 @@ git pull origin master
 
 ## How to use it?
 
-1) Optionally rename dge_master.sh to something like dge_mouse_helin.sh
-2) Adjust paths in master script
+1) Copy dge_master_template.sh into and rename into something meaningful (e.g. dge_mouse_helin.sh)
+2) Adjust DGE_HOME in master script
 3) Run on madmax
 4) Adjust reports if necessary
 

dge_workflow/dge_master.sh → dge_workflow/dge_master_template.sh

 export baseDir=<<PATH_TO_BASEDIR>>
 export project=<<PROJECTNAME>>
-
 screen -R $project
 
-########################################################################################################################
-### Setup
 
-## change if you want to use customized copy
 DGE_HOME=/projects/bioinfo/holger/bioinfo_templates/dge_workflow
-
-
-source <(curl https://dl.dropboxusercontent.com/u/113630701/datautils/bash/lsf_utils.sh 2>&1 2>/dev/null)
-source <(curl https://dl.dropboxusercontent.com/u/113630701/datautils/R/utils/spinr.sh 2>&1 2>/dev/null)
 source $DGE_HOME/dge_utils.sh
-
-
-export PATH=/projects/bioinfo/holger/bin/bowtie2-2.2.2:$PATH
-export PATH=/projects/bioinfo/holger/bin/tophat-2.0.13.Linux_x86_64:$PATH
-export PATH=/home/brandl/bin/cufflinks-2.2.1.Linux_x86_64:$PATH
-# which tophat;  which bowtie2; which cuffdiff
+export PATH=$DGE_HOME:$DGE_HOME/../misc/:$PATH
 
 
 ########################################################################################################################
-#### QC
+### QC
 
+dge_fastqc *fastq.gz &
 
-#fastqFiles=$(ls /projects/bioinfo/holger/projects/helin/all_trep_pooled/dog_*)
-fastqFiles=<<<<TBD>>>>
-#ll $fastqFiles
 
-## do some qc for the reads
-dge_fastqc -o $baseDir/fastqc $fastqFiles &
+########################################################################################################################
+### Trim low-quality bases and remove left-over adapters
+
+mcdir $baseDir/trimmed
 
-## create a fastqc report
-spinr $DGE_HOME
+dge_cutadapt $(ls $baseDir/treps_pooled/*fastq.gz) 2>&1 | tee cutadapt.log
 
-mailme "$project: fastqc done"
+dge_fastqc $(ls *fastq.gz)
 
 
 ########################################################################################################################
@@ -43,18 +29,30 @@ mailme "$project: fastqc done"
 
 mcdir $baseDir/mapped
 
-#igenome=/projects/bioinfo/igenomes/Canis_familiaris/Ensembl/CanFam3.1
+fastqFiles=$(ls $baseDir/trimmed/*.fastq.gz)
 igenome=<<<<TBD>>>>
+## Example:
+## igenome=/projects/bioinfo/igenomes/Canis_familiaris/Ensembl/CanFam3.1
 
 
 dge_tophat_se -i $igenome $fastqFiles 2>&1 | tee mapped.log
 
 mailme "$project: mapping done"
 
+
 ########################################################################################################################
-#### Basic Alginment QC and technical replicate grouping
+### Basic Alginment QC and technical replicate grouping
+
+mcdir $baseDir/trep_pooled
+
+
+bio_reps=<<<biological replicates labels>>>
+## Examples
+# bio_reps=$(csvcut  -tc bio_sample  ../lanereps_pooled/renaming_scheme.txt |  tail -n+2 | sort -u | xargs echo | tr " " ",")
+## bio_reps="ctrl,isnm1"
+
+dge_merge_treps $baseDir/mapped/ $bio_reps
 
-##TBD
 
 ########################################################################################################################
 ### Do the differential expression analysis
@@ -70,4 +68,21 @@ labels=<<<<TBD>>>>
 dge_cuffdiff $gtfFile $baseDir/mapped $labels
 MakeCuffDB $gtfFile "NAN"
 
-mailme "$project: cuffdiff done"
\ No newline at end of file
+mailme "$project: cuffdiff done"
+
+
+mcdir $baseDir/cuffdiff/dge_report
+#export DGE_PARAMS="-S"
+spin.R $DGE_HOME/dge_analysis.R $baseDir/cuffdiff
+
+
+########################################################################################################################
+### Sync back to project space
+
+# ... project specific stuff
+screen -R rsync_$project
+
+rsync -avsn --delete  $baseDir brandl@fileserver:/projects//project-sequencing-helin/rnaseq_cyst/results
+
+mailme "$project: rsync done"
+

dge_workflow/dge_utils.sh

@@ -84,6 +84,7 @@ done
 wait4jobs .cajobs
 ziprm cutadapt_logs ${project}__ca__*.log
 
+
 ## todo do a small report here about what has been trimmed away and why
 
 }
@@ -186,6 +187,53 @@ mailme "$project: bamcorrelate done in $(pwd)"
 export -f dge_bam_correlate
 
 
+## Merge technical replicatews
+dge_merge_treps(){
+
+usage="Usage: dge_bam_correlate <bam_directory> <comma_sep_biosample_spec>"
+
+if [ $# -ne 2 ]; then
+    echo $usage >&2 ; return;
+fi
+
+
+local bam_dir=$1; # bam_dir=$baseDir/mapped_trim/
+local bio_reps=$2; # bio_samples="test,lala"
+
+if [ ! -d "$bam_dir" ]; then
+    echo "bam file directory does not exist! $usage'" >&2 ; return;
+fi
+
+if [ $(echo "$bio_samples" | grep "," | wc -l )  -ne 1 ]; then
+    echo "Invalid biosample spec! $usage'" >&2 ; return;
+fi
+
+#for sample in aRG bRG N; do
+for sample in $(echo $bio_samples | tr ",", " "); do
+    # DEBUG sample=Insm1_1103
+
+    ## todo use bsub instead here
+    (
+    sampleBams=$(find $bam_dir -name '*bam'  | grep -v unmapped | grep $sample)
+    echo "merging $sample with $sampleBams"
+
+    if [ 1 -eq $(echo "$sampleBams" | wc -l) ]; then
+        cp $sampleBams $sample.bam
+    else
+        samtools merge - $(echo $sampleBams | xargs echo) | samtools sort - $sample
+    fi
+
+    samtools index $sample.bam
+    )&
+done
+
+}
+export -f dge_merge_treps
+
+## tester
+#dge_merge_treps $baseDir/mapped_trim/ "Ctrl_0703,Ctrl_1029,Ctrl_1103,Insm1_0703,Insm1_1029,Insm1_1103"
+
+
 
 dge_cuffdiff(){