-
Notifications
You must be signed in to change notification settings - Fork 0
/
bt2-map.sh
64 lines (53 loc) · 1.91 KB
/
bt2-map.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#!/bin/bash
#SBATCH --job-name=bt2map
#SBATCH --partition=amd
#SBATCH --array=1-2 #3-15,17,18,20-27,30-43,45-54 # sample numbers, currently sends only samples 1&2 to array -- modify!
#SBATCH --cpus-per-task=16
#SBATCH --mem-per-cpu=4096
#SBATCH --time=300 # Time limit hrs:min:sec
#SBATCH --output=slurm_%j.log # Standard output and error log
# exit when any command fails
set -e
# keep track of the last executed command
trap 'last_command=$current_command; current_command=$BASH_COMMAND' DEBUG
# echo an error message before exiting
trap 'echo "\"${last_command}\" command filed with exit code $?."' EXIT
eval "$(conda shell.bash hook)"
# cd to YOUR project directory, modify this line
cd /gpfs/space/home/taavi74/Projects/mapto
# slurm variables
cpus="$SLURM_CPUS_PER_TASK"
ID=EV$(printf %02d $SLURM_ARRAY_TASK_ID) # OUR sample id
# input/output variables
gb="U30316.1" # this is your ref sequence
FA=${gb}.fa
WD="results/bowtie2" # this is where your results will be
bt2base="${WD}/$(basename ${FA%.*})" #${ref%.*}
conda activate bowtie2
# QC-d reads, edit when necessary
fq1=$(ls -m /gpfs/space/projects/preterm/saja-aastased/results/trimmed_reads/${ID}*1.fastq.gz)
fq2=$(ls -m /gpfs/space/projects/preterm/saja-aastased/results/trimmed_reads/${ID}*2.fastq.gz)
input="-1 $(echo $fq1 | sed 's/ //g') -2 $(echo $fq2 | sed 's/ //g')"
# map library to reference
if [ ! -f "${bt2base}-${ID}.bam" ]
then
bowtie2 \
-p "${cpus}" \
-x "${bt2base}" \
$input \
2> "${bt2base}-${ID}.bowtie2.log" | \
samtools view -@ "${cpus}" -bS | \
samtools sort -@ "${cpus}" -o "${bt2base}-${ID}.bam"
samtools index "${bt2base}-${ID}.bam"
fi
conda deactivate
conda activate bbmap
bam="${bt2base}-${ID}.bam"
covstats="${bt2base}-${ID}.covstats.txt"
rpkm="${bt2base}-${ID}.rpkm.txt"
if [[ ! -f "$covstats" ]] | [[ ! -e "$covstats" ]]
then
echo "Working on $bam"
pileup.sh in=$bam out=$covstats rpkm=$rpkm dupecoverage=false
fi
conda deactivate