WES alignment pipeline

#!/bin/bash
# This script is to align DNA seq reads to reference genome GRCH38 and to conduct removing PCR duplicates
# The softwares and parameters in this pipeline were refered from TCGA bioinformatics pipeline
# https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/DNA_Seq_Variant_Calling_Pipeline/#dna-seq-alignment-command-line-parameters
# In case of trimming adaptor sequence, I determined to use bbmap since it will try to cut whole illumina adaptor sequence, which is possibly the most widely used adaptor sequence in the present.

# Program path
# bbmap_setting, adaptor cutting
bbmap_path=/mnt/bigHDD/resource/ref_and_tools/tools/bbmap/bbduk.sh
adaptor_seq_path=/mnt/bigHDD/resource/ref_and_tools/tools/bbmap/resources/adapters.fa

# Alinger setting
bwa_path=/mnt/bigHDD/resource/ref_and_tools/tools/bwa/bwa-0.7.17/bwa
reference_genome_path=/mnt/bigHDD/resource/ref_and_tools/reference/human/ensembl/GRCh38/Homo_sapiens.GRCh38.dna.primary_assembly.fa
reference_genome_index_path=/mnt/bigHDD/resource/ref_and_tools/reference/human/ensembl/GRCh38/Homo_sapiens.GRCh38.dna.primary_assembly.fa
gatk_dict_path=/mnt/bigHDD/resource/ref_and_tools/reference/human/ensembl/GRCh38/
# GATK setting
gatk_path=/mnt/bigHDD/resource/ref_and_tools/tools/gatk/gatk-4.0.3.0/gatk

# Samtools path
samtools_path=/mnt/bigHDD/resource/ref_and_tools/tools/samtools/samtools-1.8/samtools

# dbSNP path
dbSNP_path=/mnt/bigHDD/resource/ref_and_tools/reference/human/dbSNP/GRCh38/dbSNP150/common_all_20180418.vcf
# Read group
read_group='@RG\tID:Group1\tPL:illumina\tSM:'

# $bbmap_path in1=colo201_1.fq.gz in2=colo201_2.fq.gz out1=tmp1.fq.gz out2=tmp2.fq.gz ref=$adaptor_seq_path ktrim=r mink=11 hdist=1 tbo tpe
# Get current directory
current_path=$(pwd)
echo $current_path


#======================================================================
# Cutting adaptor sequence and alingments
#======================================================================
#===============================
# Cut adaptor sequences parameter
#===============================
echo 'Cutting adaptor sequence with bbmap program'
echo "You need to see whether bbmap's adaptor.fa has ur adaptor sequence"
mkdir 1.adaptor_cutted_fastq
outfq_dir=/1.adaptor_cutted_fastq/
outfq_path=$current_path$outfq_dir
outfq_suffix=.cutted.fq.gz

#===============================
# alingment parameters
#===============================
# Make output directory with bam files
 mkdir 2.out.DNA.bam.dir
 outbam_dir=/2.out.DNA.bam.dir/
 outbam_path=$current_path$outbam_dir
 outfile_suffix=_GRCh38_aligned.bam
 
#===========================
# Get files and pairing 1st fq and 2nd fq
#===========================
echo 'Provide file 1st suffix to recognize 1st fasta-files (ex : 1.fastq.gz)'

read suffix1

echo "The file 1st suffix is '$suffix1'"

echo 'Provide file 2nd suffix to recognize 2nd fasta-files (ex : 2.fastq.gz)'

read suffix2

echo "The file 2nd suffix is '$suffix2'"

#fasta file은 2개이므로 1P와 2P로 해둠. fasta.gz인 경우도 있지만 압축을 다 풀었다고 가정함.

# Generate input file list
ls | grep $suffix1 > input_file_1st_suffix_for_bwa.alinger.txt
ls | grep $suffix2 > input_file_2nd_suffix_for_bwa.alinger.txt
sed "s/$suffix1//" input_file_1st_suffix_for_bwa.alinger.txt > input_file_prefix_for_bwa.alinger.txt

# Generate input file pairs
paste -d',' input_file_1st_suffix_for_bwa.alinger.txt input_file_2nd_suffix_for_bwa.alinger.txt > input_file_pairs_for_bwa.alinger.txt
cat input_file_pairs_for_bwa.alinger.txt

#==========================
# If else to check whether file pair correct or not
#==========================
echo 'Is every pair correct?(y/n)'
read yes
echo "Your decided $yes"


if [ "$yes" = "y" ] || [ "$yes" = "yes" ] || [ "$yes" = "Y" ]; then

for j in $(cat input_file_pairs_for_bwa.alinger.txt)

do

echo $j

file1=$( echo $j | cut -d',' -f1)

file2=$( echo $j | cut -d',' -f2)

out_file1=$outfq_path$file1$outfq_suffix
out_file2=$outfq_path$file2$outfq_suffix
#==========================
#Do adaptor cutting process
#==========================
echo "Initiating cutting adaptor process with for loop"
$bbmap_path in1=$file1 in2=$file2 out1=$out_file1 out2=$out_file2 ref=$adaptor_seq_path ktrim=r mink=11 hdist=1 tbo tpe

#=========================================
#Do alignment process and create bam files
# Unmapped reads will be removed by samtools -F 4 function because SortSam return error for mapping quality 0 reads (MAPQ)
#=========================================
echo "Alignment reads on reference genome GRCH38"
bam_file_name_prefix=$( echo $file1 | sed "s/$suffix1//")
$bwa_path aln -t 8 $reference_genome_index_path $out_file1 > $outbam_path$(echo 'tmp1.sai')
$bwa_path aln -t 8 $reference_genome_index_path $out_file2 > $outbam_path$(echo 'tmp2.sai')
$bwa_path sampe -r $read_group$bam_file_name_prefix$(echo 'reference_GRCH38') $reference_genome_index_path $outbam_path$(echo 'tmp1.sai') $outbam_path$(echo 'tmp2.sai') $out_file1 $out_file2 | $samtools_path view -ShbF 4 -o $outbam_path$(echo 'tmp3.bam')

#=========================================
#Remove unmapped reads
#=========================================

#========================================
# Sort Bam file
#========================================
echo "Sort BAM file according to genomic coordinate"
$gatk_path SortSam --CREATE_INDEX=false \
-I=$outbam_path$(echo 'tmp3.bam')  -O=$outbam_path$(echo 'tmp4.bam') \
--SORT_ORDER=coordinate --VALIDATION_STRINGENCY=STRICT

#========================================
# Mark PCR duplicates
#========================================
echo "Mark PCR duplicates"
$gatk_path MarkDuplicates --CREATE_INDEX=true \
-I=$outbam_path$(echo 'tmp4.bam') \
-O=$outbam_path$bam_file_name_prefix$outfile_suffix \
-M=$outbam_path$(echo 'tmp4.mat')
--VALIDATION_STRINGENCY=STRICT

#========================================
# Realigntargetcreator and Indelrealigner were retired in GATK4.
# https://gatkforums.broadinstitute.org/gatk/discussion/11455/realignertargetcreator-and-indelrealigner
#========================================
# Not run realignertargetcreator-and-indelrealigner

#========================================
# BASERECALIBRATOR
# https://gatkforums.broadinstitute.org/gatk/discussion/11455/realignertargetcreator-and-indelrealigner
#========================================
echo "BASE quality recalibrator"
$gatk_path BaseRecalibrator -R $gatk_dict_path \
-I=$outbam_path$bam_file_name_prefix$outfile_suffix \
--known-sites $dbSNP_path \
-O=$outbam_path$bam_file_name_prefix$(echo '.grp')

done

fi

#========================================
# Remove temporary files
#========================================
rm -rf $outbam_path$(echo 'tmp1.sai') $outbam_path$(echo 'tmp2.sai') $outbam_path$(echo 'tmp3.bam')  $outbam_path$(echo 'tmp4.bam')

wes_alignment_pipeline.sh

0.01MB