Skip to content
Snippets Groups Projects
Commit 6fffaaee authored by elabaron's avatar elabaron
Browse files

add src/norm_coverage and add coverage calcul to RNAseq worflow

parent e6956614
No related branches found
No related tags found
No related merge requests found
......@@ -129,6 +129,7 @@ samtools index ${file_id}_sorted.bam
"""
}
sorted_bam_files.into{sorted_bam_htseq, sorted_bam_coverage}
/* HTseq */
......@@ -136,7 +137,7 @@ process sort_bam {
tag "$file_id"
input:
set file_id, file(bam) from sorted_bam_files
set file_id, file(bam) from sorted_bam_htseq
output:
set file_id, "*_htseq.bam" into sorted_bam_files_2
......@@ -188,3 +189,22 @@ htseq-count ${bam[0]} ${gtf} \
"""
}
process coverage {
tag "$file_id"
publishDir "${params.output}/05_coverage/", mode: 'copy'
input:
set file_id, file(bam) from sorted_bam_coverage
output:
file "*.bw" into coverage_files
script:
"""
bash src/norm_coverage.sh -b ${bam} \
-o {file_id}.bw \
--binSize 1 \
-p ${cpus} 8
"""
}
#!/usr/bin/env bash
set -e
usage() { echo "Usage: $0 -b <bamfile.bam> -o <outputName> --binSize <int> -p <CPUs>" 1>&2; exit 1; }
cpus=4
binSize=1
while getopts "b:o:binSize:p:" arg; do
case $arg in
-h)
echo "usage"
;;
-b)
bam=$OPTARG
;;
-o)
output=$OPTARG
;;
--binSize)
binSize=$OPTARG
;;
-p)
cpus=$OPTARG
;;
esac
done
hg38=$(samtools view ${bam} | awk '{print $1}' | sort | uniq | wc -l)
factor=$(echo "1000000/($hg38)" | bc -l)
echo "hg38 counts : $hg38"
echo "scaling factor : $factor\n"
bamCoverage -p ${cpus} --scaleFactor ${factor} --binSize ${binSize} -b ${bam} -o ${output}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment