Automation script

get_vcf.sh

#!/bin/bash
set -e
set -u
set -o pipefail


if [ "$#" -ne 1 ]; then
	echo "./get_vcf.sh "
	exit 1
fi

run_name="$1"
run_db="DRA002862.run.xml.tsv"
exp_db="DRA002862.run.xml.exp.tsv"
exp_name="$(grep $run_name $exp_db | awk -F'\t' '{ print $2 }')"
dir_name="$(grep $run_name $run_db | awk -F'\t' '{ print $3 }')"
fastq_repository="ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA002/DRA002862"

fastq1URL="$fastq_repository/$exp_name/${run_name}_1.fastq.bz2"
fastq2URL="$fastq_repository/$exp_name/${run_name}_2.fastq.bz2"

proj_dir="./stap/reads/${dir_name}"
log="get_vcf.log"

echo $dir_name

if [ ! -d "./stap" ]; then
	mkdir ./stap
fi

if [ ! -d $proj_dir ]; then
	mkdir $proj_dir
fi

cd $proj_dir

wget -c $fastq1URL
bzip2 -d ${run_name}_1.fastq.bz2
gzip ${run_name}_1.fastq
echo "got fastq file1" >> $log

wget -c $fastq2URL
bzip2 -d ${run_name}_2.fastq.bz2
gzip ${run_name}_2.fastq
echo "got fastq file2" >> $log

bwa aln -n0.02 ../../genome/GRCm38_index ${run_name}_1.fastq.gz > ${run_name}_1.sai 2>> $log
echo "got .sai file1" >> $log

bwa aln -n0.02 ../../genome/GRCm38_index ${run_name}_2.fastq.gz > ${run_name}_2.sai 2>> $log
echo "got .sai file2" >> $log

bwa sampe -a2000000 ../../genome/GRCm38_index ${run_name}_1.sai ${run_name}_2.sai ${run_name}_1.fastq.gz ${run_name}_2.fastq.gz 2>> $log | samtools view -bS - > ${run_name}.bam 2>> $log
echo "got .bam file" >> $log

samtools sort ${run_name}.bam -o ${run_name}.sort.bam 2>> $log
echo "sorted .bam file" >> $log

mv ${run_name}.sort.bam ${run_name}.bam
samtools index ${run_name}.bam
echo "indexed .bam file" >> $log

samtools mpileup -uf ../../genome/Mus_musculus.GRCm38.74.dna.toplevel.fa.gz ${run_name}.bam 2>> $log | bcftools call -vc -Oz > ${run_name}.vcf.gz 2>> $log
echo "got .vcf file" >> $log

bcftools stats ${run_name}.vcf.gz > vcfstats
echo "done" >> $log

#~/mput_hidrive.sh ${run_name}.vcf.gz
#echo "backuped .vcf" >> $log

#mkdir reports
#~/FastQC/fastqc -o reports -t 1 ${run_name}_1.fastq.gz
#~/FastQC/fastqc -o reports -t 1 ${run_name}_2.fastq.gz
#echo "generated FastQC reports" >> $log


exit 0