#fastq for all files in directory
fastqc -q *.fastq
AdapterRemoval --file1 TS-Bone-17-6_S4_L001_R1_001.fastq --basename output_single --trimns --trimqualities --gzip
bwa aln -l 1000 -n 0.01 $ref ${NAME}.trimmed.merged.fastq.gz > ${NAME}.trimmed.merged.sai   # Map with seed disabled
bwa samse $ref ${NAME}.trimmed.merged.sai ${NAME}.trimmed.merged.fastq.gz > ${NAME}.trimmed.merged.bwa.all.sam  # generate SAM alignment
  samtools view -bSh ${NAME}.trimmed.merged.bwa.all.sam > ${NAME}.trimmed.merged.bwa.all.bam  # Generate bam file
  samtools view -bh -F4 ${NAME}.trimmed.merged.bwa.all.bam > ${NAME}.trimmed.merged.mapped.bwa.bam  # Filter unmapped reads (F4 flag))
  samtools view -bh -q 30 ${NAME}.trimmed.merged.mapped.bwa.bam > ${NAME}.trimmed.merged.mapped.q30.bwa.bam  # Filter Q30 quality 
  samtools sort ${NAME}.trimmed.merged.mapped.q30.bwa.bam > ${NAME}.trimmed.merged.mapped.q30.bwa.srt.bam # Sort alignment by leftmost coordinate 
  samtools rmdup -sS ${NAME}.trimmed.merged.mapped.q30.bwa.srt.bam ${NAME}.trimmed.merged.mapped.q30.bwa.srt.rmdup.bam # Remove duplicates
  # Remove reads with multiple mappings	 
  samtools view -h ${NAME}.trimmed.merged.mapped.q30.bwa.sort.rmdup.bam | grep -v 'XT:A:R'| grep -v 'XA:Z' |grep -v 'XT:A:M' | awk '{if($0~/X1:i:0/||$0~/^@/  )print $0}' | samtools view -bS - > ${NAME}.trimmed.merged.mapped.q30.bwa.srt.rmdup.uniq.bam 


#simple variant calling
bcftools mpileup -Ou -f reference.fa ${NAME}.trimmed.merged.mapped.q30.bwa.srt.rmdup.uniq.bam | bcftools call -mv -Ob -o calls.vcf.gz
bcftools index calls.vcf.gz
bcftools norm -f reference.fa calls.vcf.gz -Ob -o calls.norm.bcf
bcftools filter --IndelGap 5 calls.norm.bcf -Ob -o calls.norm.flt-indels.bcf
cat reference.fa | bcftools consensus calls.vcf.gz > consensus.fa

samtools mpileup -uF reference.fa my.bam | bcftools call -mv -Oz -o calls.vcf.gz