Exome seq Analysis DOUBT

[kamilo889]

Hi Everyone, after surpass all the inconvenients of installation in an IMac I've finally capable to perform several steps following the GATK pipeline runing it directly in the terminal

for Illumina Data
this are the steps completed:
1. Concatenate the HG19
2. INdex the HG19 with BWA
3. Align the FASTQ sequences to the genome using BWA mem
4. Sort the .sam file to create the .bam with Picard
5. Mark duplicates with picard too
6. add the RG with picard

But I have a huge question... how I do know that my analysis is good? that my BWA are correctly aligned, that the sort and duplicates are correctly marked and the BAM file is good. exist a way to know that?
thanks a lot for your time.

Regards
Camilo Andres
Biological Engineer

[Brian Bushnell]

Try using your process on a dataset where you know the answer to confirm that it works. This could be real, publicly available data, or some synthetic data you make yourself, by for example editing a small genome (e.coli, for example) and generating synthetic reads from it, then mapping to the original genome, calling variations, and verifying that your end results detect the edit you made. If so, all the intermediate steps must have worked (aside from marking duplicates which would need to be tested differently).

[Zaag]

You could use this: http://www.bioplanet.com/gcat

[Baseless]

As above, you could use a known dataset. You could also sit down and consider what happens in which step and do some checks with samtools view from your alignment. You could just grep certain information from there, since the sam/bam format is documented. Or you could just run samtools flagstat to see how alignment went.

You might also look into introducing the RG-information at the BWA step, which saves you the I/O time to do it later with Picard. And if you are in a hurry, you might also want to sort on a datastream and not do it step by step.

The following one would take you from paired-end fastq to a sorted bam file, but be careful and watch the amount of memory consumed before you go parallel:

time bwa mem -t 8 -M -R '@RG\tID:ID_of_your_run\tCN:Name_of_your_department\tPU:ID_of_your_flowcell\tDTate_when_the_run_started\tPLevice_Platform\tLB:Library_type\tSM:SAMPLENAME' \
hGRC37.fa Sample_R1.fastq.gz Sample_R2.fastq.gz | samtools view -bSu -F 0x04 - \
| samtools sort -m 4294967296 - Sample.srt