Attention:
The frontline support team will be unavailable to answer questions until May27th 2019. We will be back soon after. Thank you for your patience and we apologize for any inconvenience!

A complete script that processes a trio with WGS data from FASTQ to BAM to VCF

Hi,

I found there are quite many GATK documentation material online. However, I could not find one that shows me exactly how to process a WGS dataset from bump to bump. For example, right now, I got 3 samples with WGS data. Each sample has 4 FASTQ files, for example for the first sample, there are: s1_L1_1.fa.gz, s1_L1_2.fa.gz, s1_L2_1.fa.gz, s1_L2_2.fa.gz. Now, my question is, what are the exact serial of commands that I should use to create a VCF file with these 3 samples.

I found a nice example at http://www.htslib.org/workflow, but I am afraid that it is not the latest version. I spent quite some time try to to figure this out. Below is what I got for running on each of the 3 samples:

  1. bwa mem -t 1000 -k 32 -M hg19.fa s1_L1_1.fa.gz s1_L1_2.fa.gz s1_L2_1.fa.gz s1_L2_2.fa.gz | samtools view -b -S -t hg19.fa.fai - > s1.bam

  2. samtools sort [email protected] 4 s1.tmp.bam s1.sorted.bam, then java -jar picard/MarkDuplicates.jar I= s1.sorted.bam O=s1.markdup.bam M=s1.dupStat

  3. gatk RealignerTargetCreator –R hg19.fq.gz –I s1.sorted.bam –known indels.vcf –O realigner.intervals

  4. gatk BaseRecalibrator –R hg19.fq.gz–I realigned.bam –knownSites dbsnp137.vcf –knownSites gold.standard.indels.vcf –O recal.table

  5. gatk HaplotypeCaller –R hg19.fa.gz –I s1.bam –o s1.gvcf –ERC GVCF

Once I done the above for each of the 3 samples, then I merge 3 gVCF files together, by:
gatk GenotypeGVCFs –R hg19.fa.gz –V s1.gvcf –V s2.gvcf –V s3.gvcf

Can someone please let me know if I got the above correct? If not, can you please kindly correct me?

Thank you & best regards,
Jie

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @jiehuang001
    Hi Jie,

    You may find this repository helpful.

    -Sheila

  • jiehuang001jiehuang001 Member

    Dear Sheila:

    Thank you very much for letting me know. I am looking for a series of BASH command, instead of WDL files. I know that WDL is supposed to be easy, but i simply want to embed all my scripts into another BASH script.

    Can you please recommend something in BASH?

    Best regards,
    Jie

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭
    edited July 2018

    Hi @jiehuang001,

    Sheila is on vacation so I am following up. You can embed BASH commands within WDL using <<<bash commands >>> heredoc syntax as described here. Cromwell and WDL allow for job management and record keeping for provenance--things researchers may be keen to outsource and automate. You can simply place your BASH commands into a WDL task.

    Note, WDL/Cromwell allows for easy parallelization across samples and genomic intervals. We also provide as well as support the provided WDL scripts (you can post questions to the WDL forum). We do not support BASH scripting. Sorry.

  • jiehuang001jiehuang001 Member

    Dear Shlee:

    Thank you very much for your reply. i will take a look at WDL then.

    Best regards,
    jie

  • jiehuang001jiehuang001 Member

    Dear Shlee:

    I still prefer to write a list of BASH commands to implement GATK, something like what i wrote in my original post. I have been used to putting all pipeline scripts in BASH commands. If you have a good reference for that, can you please share?

    Best regards,
    Jie

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @jiehuang001,

    No, we do not provide a reference for the GATK pipelines written in BASH only. However, converting the commands in the WDLs to bash should be fairly straightforward.

    Good luck!

  • jiehuang001jiehuang001 Member

    Dear all:

    Thank you for all the kind replies. In the link that Sheila originally recommended https://github.com/gatk-workflows/gatk4-germline-snps-indels, I see quite a few files. Should I first use haplotypecaller-gvcf-gatk4.wdl and then use joint-discovery-gatk4.wdl?

    What about the steps listed in the left part of the GATK best practices, including the non-GATK part ("Map to Reference" and "Mark Duplicates & Sort"), and the two steps before creating a BAM file ("Indel Realignment" and "Base Recalibration")?

    Again, I am looking for a complete series of GATK commands that process FASTQ files to VCF files, based on best practices.

    Best regards,
    Jie

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Have a look at the Best Practices pages; they include tables that point to specific WDLs for each part of the pipeline(s). https://software.broadinstitute.org/gatk/best-practices/

  • jiehuang001jiehuang001 Member

    Thanks, Geraldine! I don't see it. Can you please kindly send me the links?

    Best regards,
    Jie

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Select the pipeline you're interested in in the left hand menu then scroll down to find the table of reference implementations. Some are not yet available but those you were asking about are.

    Eg https://software.broadinstitute.org/gatk/best-practices/workflow?id=11165

  • jiehuang001jiehuang001 Member

    Hi, I still could not understand why i need to follow all these instructions on this page https://software.broadinstitute.org/wdl/documentation/article?id=7158, for example, to write at least 20 lines of code, when i simply want to call the GATK HaplotypeCaller. This should be able to be done in a single BASH line such as: gatk -T HaplotypeCaller.jar -R $ref -ERC GVCF --genotyping_mode DISCOVERY -stand_call_conf 30 -stand_emit_conf 10 -minPruning 3 -I aligned_reads.sorted.dedup.realigned.recal.BAM -o raw_variants.VCF

    Plus, i now need to install this cromwell.jar, womtool.jar, wdltool.jar, etc. Furthermore, I don't know how to embed these WDL codes into my other BASH commands and call the WDL script from other .sh scripts. So, is this WDL thing really the preferred way to run GATK? Can't I simply write a series of BASH commands to do the same thing?

    Really a lit frustrated. Can someone please help? Thanks!!!

    Jie

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @jiehuang001, I understand your frustration; let me clarify a few things.

    First, yes you can absolutely run GATK commands through BASH scripts if you prefer that over WDL. There is no requirement to use WDL, it is just a pipelining approach that we recommend because we find it convenient and highly portable. It is also what we have chosen as a way to distribute reference implementations of the GATK Best Practices, ie full pipelines, so one of the benefits of using WDL is that if you just want our standard pipelines you can just use what we provide and you don’t have to write it all yourself (or you can use them as a base and modify whatever you want to do differently).

    To be clear, we would not expect you to embed WDL within your existing BASH scripts — WDL + Cromwell are intended to replace those entirely.

    Certainly if you’re only running one command at a time the WDL form seems like overkill over a simpler invocation through BASH, but when you’re running full pipelines, WDL and Cromwell give you a lot of flexibility for doing things like complex flow control, smart resume etc that are much harder to do in pure BASH.

    So ultimately it’s up to you to choose how you prefer to work, and how much you’re willing to do yourself vs taking an off the shelf solution.

Sign In or Register to comment.