A bash script for running on a bunch of bam files

Question

I have some bam files in this directory

/data/Continuum/WES/results/

I want to run GATK mutation calling over bam files

I googled and I realised for one function, I can do this

cd /data/Continuum/WES/vcf/

for file in *.bam ; do ./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta done;

I have five functions to run but I do not know how to include all in the script

These are my functions

./gatk CollectSequencingArtifactMetrics -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta


./gatk GetPileupSummaries -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -V af-only-gnomad.hg38.vcf.gz -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta


./gatk CalculateContamination -I NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.table -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.contamination.table

./gatk GetSampleName -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name

./gatk Mutect2 -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta -I /data/Continuum/WES/results/NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.bam -O NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.vcf -tumor NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup.targeted_sequencing.sample_name --af-of-alleles-not-in-resource 2.5e-06 --germline-resource af-only-gnomad.hg38.vcf.gz -pon 1000g_pon.hg38.vcf.gz

I thought to use | to include five functions but I do not know how to change input and output for each iteration

I tried this for the first function,

    FILES="/data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*"
for f in $FILES
do
  ./gatk CollectSequencingArtifactMetrics -I $f -O $f --FILE_EXTENSION .txt -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta
done

but returned this error

SAMException: Cannot read non-existent file: file:///data/Continuum/WES/testAlignmentBROADGenome2/results/.MarkedDup.bam*

Anyone, please help me in develop this?

Answer 1

There are a couple ways to approach this, one is to create output files based on the BAM input file name - the other approach is to make a new folder (based on the input BAM filename containing result files with the same name for each input. I generally use the first method, but either work.

Here is a script that may run the 5 GATK functions. It calculates a prefix to be used for derivative output files BASE, and should run your pipeline on all the BAM files within the INPUT_FOLDER.

#!/bin/bash

set -e

INPUT_FOLDER="/data/Continuum/WES/results/"

for BAM_FILE in "$INPUT_FOLDER/"*.bam
do
        BASE=$(basename "$BAM_FILE" .bam) # NG-27280_CLTSS_LTS_0017_lib506243_7661_2_MarkedDup

        # function 1
        ./gatk CollectSequencingArtifactMetrics \
                -I "$BAM_FILE" \
                -O "$BASE" \
                --FILE_EXTENSION .txt \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 2
        ./gatk GetPileupSummaries \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.table" \
                -V af-only-gnomad.hg38.vcf.gz \
                -L resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta

        # function 3
        ./gatk CalculateContamination \
                -I "$BASE.targeted_sequencing.table" \
                -O "$BASE.targeted_sequencing.contamination.table" \

        # function 4
        ./gatk GetSampleName \
                -I "$BAM_FILE" \
                -O "$BASE.targeted_sequencing.sample_name"

        # function 5
        ./gatk Mutect2 \
                -R resources_broad_hg38_v0_Homo_sapiens_assembly38.fasta \
                -I "$BAM_FILE" \
                -O "$BASE".vcf \
                -tumor "$BASE.targeted_sequencing.sample_name" \
                --af-of-alleles-not-in-resource 2.5e-06 \
                --germline-resource af-only-gnomad.hg38.vcf.gz \
                -pon 1000g_pon.hg38.vcf.gz
done

It may need some adjustments - I did not test this, and used the GATK incantations you provided (thank you for including all of your steps!)

It includes the addition of the command set -e, which depending on GATKs return status may stop the script if any command fails - otherwise the script will happily run through a failed command.

You also may need to adjust the location of your output, and the directory the script is in upon execution, as you are running GATK from the current directory (./gatk)

In order to make these long lines readable, I use the method of placing a \ at the end of a line - this allows the command to wrap to the next line. Be sure that the \ character is at the end of your line for this to function properly.