-
Notifications
You must be signed in to change notification settings - Fork 0
/
workflow1.sh
53 lines (37 loc) · 1.35 KB
/
workflow1.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
#!/bin/bash
# Updates to work with samtools Version: 0.1.16 (r963:234)
# and bwa Version: 0.6.1-r104
# and bedtools v2.16.2
#input is a fastq file '$1'
export DB=/homes/genomes/s.cerevisiae/sacCer3/bwa_indexes/sacCer3
#bwa to align
echo "## Aligning using bwa"
bwa aln $DB $1 -t 40 > $1.sai
bwa samse $DB $1.sai $1 > $2.sam
#sam to bam
echo "## Converting sam to bam, and remove unmapped reads"
samtools view -h -b -S -F 4 $2.sam > $2.bam
#sort
echo "## sorting bam"
samtools sort -m 500000000000 $2.bam $2-all
#remove 'dupicates'
echo "## removing duplicates"
samtools rmdup -s $2-all.bam $2-rmdup.bam
#remove intermediate #tidy up
echo "## remove intermediate & tidy up"
rm $2.bam
gzip -9 $2.sam &
#remove rDNA region
echo "## intersecting"
bedtools intersect -v -abam $2-rmdup.bam -b rDNA.bed > $2-rDNA-rmdup.bam
bedtools intersect -v -abam $2-all.bam -b rDNA.bed > $2-rDNA-all.bam
#Info for rDNA alone
#this was ignored as region nr saturated with reads
echo "## generating coverage bed"
bedtools coverage -abam $2-rmdup.bam -b rDNA.bed -d > $2-rmdup-rDNA.bed
#This used in figures
bedtools coverage -abam $2-all.bam -b rDNA.bed -d > $2-all-rDNA.bed
### Additional steps
## To create bigwig in RPM for IGV:
# genomeCoverageBed -ibam "bamfile" -bg -scale "1000000/mapped reads" > file.bg
# wigToBigWig "file.bg" sacCer3.fa.lengths file.bigwig