0
$\begingroup$

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?

$\endgroup$

1 Answer 1

6
$\begingroup$

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.

$\endgroup$
6
  • 1
    $\begingroup$ Wow, cool, thanks. $\endgroup$
    – M__
    Commented Aug 21, 2023 at 3:09
  • $\begingroup$ OMG Thanks a million the script is running $\endgroup$
    – Zizogolu
    Commented Aug 21, 2023 at 10:21
  • $\begingroup$ Sorry @Scot, as the Mutect2 part of the script takes 15 hours in average for each sample, does this mean that for function 5, script will wait for the first sample to being completed by Mutect2 then goes for the second sample or Mutect2 will run of all of 18 samples in parallel? In the first case, the script will take 18 days? $\endgroup$
    – Zizogolu
    Commented Aug 21, 2023 at 10:27
  • $\begingroup$ Shell scripts are run line by line. If a line takes 15 hours, the interpreter will not move to the next line before that line is done. Loops are merely instructions to run a set of lines multiple times in series, so there's no parallel action happening here. You should ideally write a snakemake script with a rule for each GATK step then run a bunch of snakemake processes in parallel using whatever parallelization means are available to you. $\endgroup$
    – Ram RS
    Commented Aug 21, 2023 at 15:24
  • 1
    $\begingroup$ Hi @Angel - yes, the current version of this script runs sequentially, and is not parallelized. If function5 Mutect2 takes a long time, then it may be worthwhile to pull it out of the loop, and run it upon completion of all other steps, and addition, explore if it can be run in parallel. $\endgroup$
    – Scot
    Commented Aug 21, 2023 at 17:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.