Genome Analysis Toolkit

Ramesh · August 26, 2024 16:16

Hello ,

I very recently started my journey in Bioinformatics.

I have been tasked with building a variant calling pipeline (WES & WGS). So far I have successfully set it up for calling SNPs and INDELs using GATKs best practices workflow in normal and dragen mode. I am trying to proceed with CNVs for whole genome (hg19). I am not able to find the right interval list from GATK's public resource bundle. I used the command:

gatk PreprocessIntervals         -R Homo_sapiens_assembly19.fasta         --padding 0         -imr OVERLAPPING_ONLY         -O hg19.preprocessed.interval_list

to generate my own interval list. I proceeded with the steps in https://gatk.broadinstitute.org/hc/en-us/articles/360035531152--How-to-Call-rare-germline-copy-number-variants#1 and built the below sequence of commands and tested it on a batch of 3 samples (I've added only the rough version without path or variable information here). I am not sure if I am doing it right. I build a cohort model from my 3 samples and called each of the samples again in case mode against the same model. Where do I use the panel of normals I created? And since I am doing this for the entire genome without specific intervals, do I skip the

--output-genotyped-segments

flag in the Germline CNV caller. Also how do I filter the CNVs after calling them? Is there a standard filtration protocol like we do for SNPs and INDELs?"

# ------------------------------------------------------------------------------------------------------------------

# STEP 1: Collect Read Counts & Create Panel of Normals- CNV

# ------------------------------------------------------------------------------------------------------------------

 $gatk PreprocessIntervals \

        -R ${ref} \

        --padding 0 \

        -imr OVERLAPPING_ONLY \

        -O ${interval}




# Collect read counts for each BAM file

cd ${aligned_reads}



# Find all bam files in given Sample directory

files=$(find . -type f -name "*GATK_sorted_dedup_reads.bam")

for i in ${files}; do

    filename=$(basename "$i" .bam)  # Extract the base filename without the .bam extension

    $gatk CollectReadCounts \

        -I $i \

        --interval-merging-rule OVERLAPPING_ONLY \

        -L ${interval} \

        -O "${data}/${filename}_read_counts.hdf5"



$gatk AnnotateIntervals \

        -L ${interval} \

        -R ${ref} \

        -imr OVERLAPPING_ONLY \

        -O ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_annotated_interval_list.tsv



SAMPLE_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    SAMPLE_INPUT+="--input $counts_file "

done




$gatk FilterIntervals \

        -L ${interval} \

        --annotated-intervals ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_annotated_interval_list.tsv \

        $SAMPLE_INPUT \

        -imr OVERLAPPING_ONLY \

        -O ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_filtered.interval_list



SAMPLE_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    SAMPLE_INPUT+="--input $counts_file "

done

$gatk CreateReadCountPanelOfNormals \

  $SAMPLE_INPUT \

-O pon.hdf5



# ------------------------------------------------------------------------------------------------------------------

# STEP 2: Denoise Read Counts

# ------------------------------------------------------------------------------------------------------------------



for counts_file in ${data}/*_read_counts.hdf5; do

    filename=$(basename "$counts_file" "_read_counts.hdf5")  # Extract the base filename without the read_counts.hdf5 extension

    $gatk DenoiseReadCounts \

        -I $counts_file \

        --count-panel-of-normals pon.hdf5 \

        --standardized-copy-ratios "${data}/${filename}_standardized_CR.tsv" \

        --denoised-copy-ratios "${data}/${filename}_denoised_CR.tsv"




    $gatk PlotDenoisedCopyRatios \

          --standardized-copy-ratios "${data}/${filename}_standardized_CR.tsv" \

          --denoised-copy-ratios "${data}/${filename}_denoised_CR.tsv" \

          --sequence-dictionary ${dict} \

          --output-prefix ${filename} \

          -O .

 
done
# ------------------------------------------------------------------------------------------------------------------

# STEP 3: Determine Germline Contig Ploidy

# ------------------------------------------------------------------------------------------------------------------



PLOIDY_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    PLOIDY_INPUT+="--input $counts_file "

done



$gatk DetermineGermlineContigPloidy \

    -L ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_filtered.interval_list \

    --interval-merging-rule OVERLAPPING_ONLY \

    $PLOIDY_INPUT \

    --contig-ploidy-priors "${ploidy_priors}" \

    --output ${data} \

    --output-prefix ploidy \

    --verbosity DEBUG



CPLOIDY_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    CPLOIDY_INPUT+="--input $counts_file "

done




$gatk DetermineGermlineContigPloidy \

    --model ${data}/ploidy-model \

    $CPLOIDY_INPUT \

    --output ${data} \

    --output-prefix ploidy_case \

    --verbosity DEBUG



# ------------------------------------------------------------------------------------------------------------------

# STEP 3: Run GermlineCNVCaller in Cohort Mode

# ------------------------------------------------------------------------------------------------------------------



COHORT_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    COHORT_INPUT+="--input $counts_file "

done

$gatk GermlineCNVCaller \

    -L ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_filtered.interval_list \

    --interval-merging-rule OVERLAPPING_ONLY \

    --run-mode COHORT \

    --contig-ploidy-calls ${data}/ploidy-calls \

    $COHORT_INPUT \

    --annotated-intervals ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19_annotated_interval_list.tsv \

    --output ${data}/CNVcalls \

    --output-prefix cohort \

    --verbosity DEBUG

# ------------------------------------------------------------------------------------------------------------------

# STEP 4: Run GermlineCNVCaller in Case Mode for Each Sample

# ------------------------------------------------------------------------------------------------------------------



CASE_INPUT=""

for counts_file in ${data}/*_read_counts.hdf5; do

    CASE_INPUT+="--input $counts_file "

done




#    filename=$(basename "$counts_file" _read_counts.hdf5)

    $gatk GermlineCNVCaller \

        --run-mode CASE \

        --contig-ploidy-calls ${data}/ploidy_case-calls \

        --model ${data}/CNVcalls/cohort-model \

        $CASE_INPUT \

        --output ${data}/CNVcasecalls  \

        --output-prefix pcase \

        --verbosity DEBUG

   

# ------------------------------------------------------------------------------------------------------------------

# STEP 5: Postprocess Germline CNV Calls to Generate VCF Files

# ------------------------------------------------------------------------------------------------------------------



i=0

for counts_file in ${data}/*_read_counts.hdf5; do

    filename=$(basename "$counts_file" _read_counts.hdf5)

    echo "Iteration $i"

    echo "Generating final vcf file for: $counts_file"

    $gatk PostprocessGermlineCNVCalls \

        --model-shard-path ${data}/CNVcalls/cohort-model \

        --calls-shard-path ${data}/CNVcalls/cohort-calls \

        --contig-ploidy-calls ${data}/ploidy-calls \

        --sample-index $i\

        --sequence-dictionary ${MAIN_DIR}/SUPP_FILES/REFERENCES/${reference}/Homo_sapiens_assembly19.dict \

        --output-genotyped-intervals ${results}/${filename}_gen_intv.vcf.gz \

        --output-genotyped-segments ${results}/${filename}_gen_seg.vcf.gz \

        --output-denoised-copy-ratios ${results}/${filename}_denoised_copy_ratios.tsv

 

    i=$((i+1))

done

Thank you, I really appreciate your help in this!

Genome Analysis Toolkit

Need Help?

Community Forum

GATK Germline CNV calling

8 comments

Welcome

Didn't find what you were looking for?

Quick Links

Recent GATK News

About the GATK community