-
Notifications
You must be signed in to change notification settings - Fork 3
/
bwa_pipeline.sh
executable file
·112 lines (82 loc) · 4.88 KB
/
bwa_pipeline.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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#!/usr/bin/bash
######### Local paths
WD=/N/u/rtraborn/Carbonate/scratch/DaphniaVariantCall
######### Loading carbonate-specific resources (change if not using this cluster) #######
module load java
#module load trimmomatic #this loads version 0.36
module load fastqc
module load bwa
######### Paths to binaries #########
samtools=/N/soft/rhel6/samtools/1.3.1/bin/samtools
picard='java -jar /N/soft/rhel6/picard/2.8.1/picard.jar'
Trimmomatic='java -jar /N/dc2/projects/daphpops/Software/Trimmomatic-0.36/trimmomatic-0.36.jar'
GATK='java -jar /N/soft/rhel6/gatk/3.4-0/GenomeAnalysisTK.jar'
bamUtil=/N/soft/rhel6/bamUtil/1.0.13/bam
# please install ngsutils (see INSTALL.txt for instructions)
# then, provide the path to the fastqutils binary
# fastqutils=/N/u/rtraborn/Carbonate/scratch/DaphniaVariantCall/software/ngsutils/bin/fastqutils
####### Path to adapters #######
adapterTrim=../adapters/Bioo_Adapters.fa
###### Misc names ######
assemblyID=PA42_4
assemblyName=PA42_with_mt.fasta
sequenceDict=PA42_with_mt.dict
##### Sequence names #####
clone1_R1=KAP-00074_CGCTGATC_L008_R2_001.fastq
clone1_R2=KAP-00074_CGCTGATC_L008_R2_001.fastq
cd $WD
echo "Copying sample files (from clone KAP-00074) to a new directory called fastq/."
mkdir fastq
cp /N/dc2/projects/daphpops/Population_samples/KAP2013/140501/Sample_KAP-00074/*.fastq fastq
echo "Creating a symbolic link to the PA42 assembly"
mkdir assembly
cd assembly
ln -s /N/dc2/projects/daphpops/PA42_with_mt/PA42_with_mt.fasta $assemblyName
echo "Making a samtools index file"
$samtools faidx $assemblyName
echo "Making a dictionary file of a reference."
# The Java version should not matter as long as it works. (ed: ?)
$picard CreateSequenceDictionary R=$assemblyName O=$sequenceDict
# 0. Performing fastqc on sample
echo "Performing fastqc on read pairs for this sample."
cd ../fastq
fastqc $clone_R1 $clone_R2
# 1. After preparing the FASTA file of adapter sequences, trim adapter sequences from sequence reads.
echo "Trimming adapter sequences from sequence reads."
$Trimmomatic PE $clone1_R1 $clone1_R2 KAP-00074_15lanes_R1-paired.fastq KAP-00074_15lanes_R1-unpaired.fastq KAP-00074_15lanes_R2-paired.fastq KAP-00074_15lanes_R2-unpaired.fastq HEADCROP:3 ILLUMINACLIP:$adapterTrim:2:30:10:2 SLIDINGWINDOW:4:15 MINLEN:30
# The additional information on Trimmomatic arguments was omitted because it is available in the Trimmomatic documentation.
# 2. Building the Hisat2 index using the assembly
bwa index ../assembly/$assemblyName
# 3. Map reads to the reference sequence using Hisat2.
echo "Mapping reads to the reference genome using bwa."
bwa mem -t 8 ../assembly/$assemblyName KAP-00074_15lanes_R1-paired.fastq KAP-00074_15lanes_R2-paired.fastq > KAP-00074_PA42_with_mt.sam
# 4. Convert the SAM file to the BAM file.
echo "Converting the sam file to bam."
$samtools view -bS KAP-00074_PA42_with_mt.sam > KAP-00074_PA42_with_mt.bam
# 5. Sort the BAM file using Picard.
echo "Sorting the bam file using Picard."
$picard SortSam INPUT=KAP-00074_PA42_with_mt.bam OUTPUT=Sorted_KAP-00074_PA42_with_mt.bam SORT_ORDER=coordinate
# 6. Add read groups to the sorted BAM file.
echo "Adding read groups to the sorted bam file."
$picard AddOrReplaceReadGroups INPUT=Sorted_KAP-00074_PA42_with_mt.bam OUTPUT=RG_Sorted_KAP-00074_PA42_with_mt.bam RGID=Daphnia RGLB=bar RGPL=illumina RGSM=KAP-00074 RGPU=6
# 7. Mark duplicate reads.
echo "Marking duplicates using Picard."
$picard MarkDuplicates INPUT=RG_Sorted_KAP-00074_PA42_with_mt.bam OUTPUT=dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam METRICS_FILE=KAP-00074_PA42_with_mt_metrics.txt
# 8. Index the BAM file using Picard.
echo "Indexing the bam file using Picard."
$picard BuildBamIndex INPUT=dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam
# 9. Define intervals to target for the local realignment.
echo "Defining intervals to target for local realignment using Picard."
$GATK -T RealignerTargetCreator -R ../assembly/$assemblyName -I dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam -o KAP-00074_PA42_with_mt.intervals
# 10. Locally realign reads around indels.
echo "Performing the local realignment using Picard."
$GATK -T IndelRealigner -R ../assembly/$assemblyName -I dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam -targetIntervals KAP-00074_PA42_with_mt.intervals -o realigned_dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam
# 11. Clip overlapping read pairs.
echo "Clipping the overlapping read pairs using bamUtil."
$bamUtil clipOverlap --in realigned_dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam --out Clipped_realigned_dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam
# 12. Index the clipped BAM file using Samtools
echo "Indexing the clipped BAM file using Samtools."
$samtools index Clipped_realigned_dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam
# 13. Make the mpileup file from the BAM file.
echo "Creating the mpileup file from the BAM file."
$samtools mpileup -f ../assembly/$assemblyName Clipped_realigned_dedup_RG_Sorted_KAP-00074_PA42_with_mt.bam -o KAP-00074_PA42_with_mt.mpileup