DNA Alignment

Introduction

Tutorials

Command Line

Setup

Source Code and Test Data

cd ~/git

git clone [email protected]:genome/analysis-workflows.git

cd analysis-workflows/

Docker Image

docker run -v `pwd`:/analysis-workflows -it mgibio/dna-alignment /bin/bash -l

Unpack Test Data Bundle

cd /analysis-workflows/example_data/

tar -xzvf exome_workflow.tgz

cd exome_workflow

Inputs

# How were these inputs created?
cat README

Reference

wc chr17_test.fa
 
head chr17_test.fa

grep -v N chr17_test.fa | head

grep -v ">" chr17_test.fa | wc

cat chr17_test.fa | grep -v ">" | perl -ne 'chomp $_; $bases{$_}++ for split //; if (eof){print "$_ $bases{$_}\n" for sort keys %bases}'

BWA Index

OPTIONAL: Takes 5 minutes

rm chr17_test.fa.*
/usr/local/bin/bwa index chr17_test.fa

SAMTools Index

rm chr17_test.fa.fai 
/opt/samtools/bin/samtools faidx chr17_test.fa
head chr17_test.fa.fai

Sequence Dictionary

rm chr17_test.dict
/usr/bin/java -jar /opt/picard/picard.jar CreateSequenceDictionary R=chr17_test.fa O=chr17_test.dict
cat chr17_test.dict

Known Variation

ls -1 *.vcf.gz
zcat chr17_test_dbsnp.vcf.gz | grep -v '^##' | head

Interval List

ls -1 *.interval_list

cat chr17_test_genes.interval_list

Unaligned BAM

/opt/samtools/bin/samtools view -H 2895499399.bam
/opt/samtools/bin/samtools view 2895499399.bam | head

Convert Unaligned BAM to FASTQ

/usr/bin/java -jar /opt/picard/picard.jar SamToFastq I=2895499399.bam INCLUDE_NON_PF_READS=true F=2895499399_1.fastq.gz F2=2895499399_2.fastq.gz

zcat 2895499399_1.fastq.gz | head
zcat 2895499399_2.fastq.gz | head

/usr/bin/java -jar /opt/picard/picard.jar SamToFastq I=2895499331.bam INCLUDE_NON_PF_READS=true F=2895499331_1.fastq.gz F2=2895499331_2.fastq.gz

Alignment

BWA MEM

/usr/local/bin/bwa mem -K 100000000 -t 2 -Y -R "@RG\tID:2895499399\tPU:H7HY2CCXX.4.TGACCACG\tSM:H_NJ-HCC1395-HCC1395_BL\tLB:H_NJ-HCC1395-HCC1395_BL-lg21-lib1\tPL:Illumina\tCN:WUGSC" chr17_test.fa 2895499399_1.fastq.gz 2895499399_2.fastq.gz > 2895499399.sam

/usr/local/bin/bwa mem -K 100000000 -t 2 -Y -R "@RG\tID:2895499331\tPU:H7HY2CCXX.3.TGACCACG\tSM:H_NJ-HCC1395-HCC1395_BL\tLB:H_NJ-HCC1395-HCC1395_BL-lg21-lib1\tPL:Illumina\tCN:WUGSC" chr17_test.fa 2895499331_1.fastq.gz 2895499331_2.fastq.gz > 2895499331.sam

addMate Tags

Add mate cigar and mate quality tags

/usr/local/bin/samblaster -a --addMateTags -i 2895499399.sam -o 2895499399_mateTag.sam
/usr/local/bin/samblaster -a --addMateTags -i 2895499331.sam -o 2895499331_mateTag.sam

SAM to BAM

/opt/samtools/bin/samtools view -b -S -o 2895499399_mateTag.bam 2895499399_mateTag.sam
/opt/samtools/bin/samtools view -b -S -o 2895499331_mateTag.bam 2895499331_mateTag.sam

alignment_helper.sh

Align the normal sample again with the helper script. NOTE: This is the same process as previous steps but using unaligned BAM files directly as the input. See recommendation #10 from this article: https://jmd.amjpathol.org/article/S1525-1578(17)30373-2/fulltext

/bin/bash /usr/bin/alignment_helper.sh 2895499399.bam "@RG\tID:2895499399\tPU:H7HY2CCXX.4.TGACCACG\tSM:H_NJ-HCC1395-HCC1395_BL\tLB:H_NJ-HCC1395-HCC1395_BL-lg21-lib1\tPL:Illumina\tCN:WUGSC" chr17_test.fa > 2895499399_alignment_helper.bam

/bin/bash /usr/bin/alignment_helper.sh 2895499331.bam "@RG\tID:2895499331\tPU:H7HY2CCXX.3.TGACCACG\tSM:H_NJ-HCC1395-HCC1395_BL\tLB:H_NJ-HCC1395-HCC1395_BL-lg21-lib1\tPL:Illumina\tCN:WUGSC " chr17_test.fa > 2895499331_alignment_helper.bam

Align a different sample, ex. tumor, using the helper script.

/bin/bash /usr/bin/alignment_helper.sh 2895499237.bam "@RG\tID:2895499237\tPU:H7HY2CCXX.4.ATCACGGT\tSM:H_NJ-HCC1395-HCC1395\tLB:H_NJ-HCC1395-HCC1395-lg24-lib1\tPL:Illumina\tCN:WUGSC" chr17_test.fa > 2895499237_alignment_helper.bam

/bin/bash /usr/bin/alignment_helper.sh 2895499223.bam "@RG\tID:2895499223\tPU:H7HY2CCXX.3.ATCACGGT\tSM:H_NJ-HCC1395-HCC1395\tLB:H_NJ-HCC1395-HCC1395-lg24-lib1\tPL:Illumina\tCN:WUGSC" chr17_test.fa > 2895499223_alignment_helper.bam

Merge

Merge Read Group BAMs

/opt/samtools/bin/samtools merge H_NJ-HCC1395-HCC1395_BL_merged.bam 2895499399_alignment_helper.bam 2895499331_alignment_helper.bam

/opt/samtools/bin/samtools merge H_NJ-HCC1395-HCC1395_merged.bam 2895499237_alignment_helper.bam 2895499223_alignment_helper.bam

Duplicate Marking

Query Name Sort

/usr/bin/sambamba sort -t 2 -n -o H_NJ-HCC1395-HCC1395_BL_nameSorted.bam H_NJ-HCC1395-HCC1395_BL_merged.bam

/usr/bin/sambamba sort -t 2 -n -o H_NJ-HCC1395-HCC1395_nameSorted.bam H_NJ-HCC1395-HCC1395_merged.bam

Mark Duplicates

/usr/bin/java -jar /opt/picard/picard.jar MarkDuplicates I=H_NJ-HCC1395-HCC1395_BL_nameSorted.bam O=H_NJ-HCC1395-HCC1395_BL_markDups.bam ASSUME_SORT_ORDER=queryname METRICS_FILE=H_NJ-HCC1395-HCC1395_BL_mark_dups_metrics.txt QUIET=true COMPRESSION_LEVEL=0 VALIDATION_STRINGENCY=LENIENT

/usr/bin/java -jar /opt/picard/picard.jar MarkDuplicates I=H_NJ-HCC1395-HCC1395_nameSorted.bam O=H_NJ-HCC1395-HCC1395_markDups.bam ASSUME_SORT_ORDER=queryname METRICS_FILE=H_NJ-HCC1395-HCC1395_mark_dups_metrics.txt QUIET=true COMPRESSION_LEVEL=0 VALIDATION_STRINGENCY=LENIENT

Coordinate Sort

/usr/bin/sambamba sort -t 2 -o H_NJ-HCC1395-HCC1395_BL_coordSorted.bam H_NJ-HCC1395-HCC1395_BL_markDups.bam

/usr/bin/sambamba sort -t 2 -o H_NJ-HCC1395-HCC1395_coordSorted.bam H_NJ-HCC1395-HCC1395_markDups.bam

Base Quality Recalibration

Calculate BQSR

/usr/bin/java -jar /opt/GenomeAnalysisTK.jar -T BaseRecalibrator -o H_NJ-HCC1395-HCC1395_BL_bqsr.table --preserve_qscores_less_than 6 --disable_auto_index_creation_and_locking_when_reading_rods --disable_bam_indexing -dfrac .1 -nct 2 -R chr17_test.fa -I H_NJ-HCC1395-HCC1395_BL_coordSorted.bam -L chr17 -knownSites chr17_test_dbsnp.vcf.gz -knownSites chr17_test_known_indels.vcf.gz -knownSites chr17_test_mills.vcf.gz

/usr/bin/java -jar /opt/GenomeAnalysisTK.jar -T BaseRecalibrator -o H_NJ-HCC1395-HCC1395_bqsr.table --preserve_qscores_less_than 6 --disable_auto_index_creation_and_locking_when_reading_rods --disable_bam_indexing -dfrac .1 -nct 2 -R chr17_test.fa -I H_NJ-HCC1395-HCC1395_coordSorted.bam -L chr17 -knownSites chr17_test_dbsnp.vcf.gz -knownSites chr17_test_known_indels.vcf.gz -knownSites chr17_test_mills.vcf.gz

Apply BQSR

/usr/bin/java -jar /opt/GenomeAnalysisTK.jar -T PrintReads -o H_NJ-HCC1395-HCC1395_BL_bqsr.bam -preserveQ 6 -SQQ 10 -SQQ 20 -SQQ 30 -nct 2 --disable_indel_quals -R chr17_test.fa -I H_NJ-HCC1395-HCC1395_BL_coordSorted.bam -BQSR H_NJ-HCC1395-HCC1395_BL_bqsr.table

/usr/bin/java -jar /opt/GenomeAnalysisTK.jar -T PrintReads -o H_NJ-HCC1395-HCC1395_bqsr.bam -preserveQ 6 -SQQ 10 -SQQ 20 -SQQ 30 -nct 2 --disable_indel_quals -R chr17_test.fa -I H_NJ-HCC1395-HCC1395_coordSorted.bam -BQSR H_NJ-HCC1395-HCC1395_bqsr.table

CRAM

Convert to CRAM

/opt/samtools/bin/samtools view -C -T chr17_test.fa -o H_NJ-HCC1395-HCC1395_BL.cram H_NJ-HCC1395-HCC1395_BL_bqsr.bam

/opt/samtools/bin/samtools view -C -T chr17_test.fa -o H_NJ-HCC1395-HCC1395.cram H_NJ-HCC1395-HCC1395_bqsr.bam

Index CRAM

/opt/samtools/bin/samtools index H_NJ-HCC1395-HCC1395_BL.cram

/opt/samtools/bin/samtools index H_NJ-HCC1395-HCC1395.cram

For additional details and a WGS tutorial: Precision Medicine Workshop: Alignment

CWL Workflow

Steps

INSERT PROCESS DIAGRAM

Inputs

Name	Description	Example	Required
bqsr_intervals
bams			✓
dbsnp_vcf			✓
known_indels			✓
mills			✓
readgroups			✓
reference			✓

DNA Alignment

Introduction

Tutorials

Command Line

Setup

Source Code and Test Data

Docker Image

Unpack Test Data Bundle

Inputs

Reference

BWA Index

SAMTools Index

Sequence Dictionary

Known Variation

Interval List

Unaligned BAM

Convert Unaligned BAM to FASTQ

Alignment

BWA MEM

addMate Tags

SAM to BAM

alignment_helper.sh

Merge

Merge Read Group BAMs

Duplicate Marking

Query Name Sort

Mark Duplicates

Coordinate Sort

Base Quality Recalibration

Calculate BQSR

Apply BQSR

CRAM

Convert to CRAM

Index CRAM

CWL Workflow

Steps

Inputs

Outputs

Known Issues

Want to contribute to this Wiki?

Clone this wiki locally