This readme describes the alignment data available on the ftp site, how it is
processed and what summary information is available for each alignment.
A. File format and alignment index
All alignment data is in the BAM format. This is the binary version of the SAM
format which is described here: http://samtools.sourceforge.net/. Each BAM file
has two associated files: a index file which has the same name but ends with
.bai and a statistics file which has the same name but end with .bas. The format
of bas file is described further down in in section D.
The alignments are found under data/XXXXXXXX/alignment and data/XXXXXXX/exome_alignment
where the XXXXXXXX identifier is the sample name.
There is an alignment.index and an exome.alignment.index which links the alignment bam
files together with their matching bai file and bas file and gives md5sums for each.
The columns are as follows:
1. bam file
2. MD5
3. bai index file
4. bai index file md5
5. bas statistics file
6. bas statistic file md5
Directory alignment_indices/ contains the most current alignment.index (the one
with the latest date and identical to the alignment.index one level up), as well
as alignment.index files of previous BAM releases. The directory also contains
the following files about BAM statistics:
yyyymmdd.alignment.index.bas.gz - a collective bas file of all BAM files.
yyyymmdd_yyyymmdd.alignment_stats.csv - a summary statistics of BAMs in any
two releases as specified by the dates in the file name; a comparison of
the two releases is captured in the "diff" values of the file. The file
contains the following information break down by platform.
1. mapped basepairs in Gb
2. # of individuals in the release
3. # of individuals with > 10 Gb of mapped sequences
mapped basepairs in Gb is also shown as break down by population at the
bottom of the file.
20091216.alignment_stats.csv - a summary statistics of the very first
release of main project BAMs.
The exome alignments also have a HsMetrics files which contain the results from
the picard tool CalculateHsMetrics
The bam filenames themselves contain a lot of information, e.g:
NA12878.mapped.ILLUMINA.bwa.CEU.low_coverage.20121211.bam
The name can be broken down into 7 pieces:
1. Sample name: this matches column 10 in the sequence.index file and
represents the individual all the sequences belong to.
2. Region, this is generally either mapped or unmapped, the mapped files represents
all the mapping to the reference genome, the unmapped file represents all the unmapped reads.
You may also see chromN files which represent mappings to just that chromosome.
3. Sequencing platform: all our current release should be ILLUMINA but old alignment files
will have SOLID or LS454.
4. Mapping algorithm: For ILLUMINA this is bwa, again older files may have bfast for solid
or ssaha or smalt for 454.
5. Population*: the Sample population 3 letter code, this code is defined
in README.populations
6. Analysis group*: this describes the sequencing strategy used for the
sequence being aligned, those strategies being 'low coverage',
'high coverage', 'exon targeted' and 'exome'.
7. Date in the format YYYYMMDD: this should match the date from the
sequence.index file which was used to produce the alignments. This date
will also be in the alignment index name. Over time one release of alignments
will contain multiple index dates. The index itself is named for the most recent
sequence.index it is based on. The index in the bam file name reflects the last
sequence index new data was added for that individual
*Note: For historical reason, in the first release of BAM files, fields 5 and 6
are replaced with a single column of project name such as SRP000033. One example
is NA20828.chrom7.ILLUMINA.bwa.SRP000033.20091216.bam
B. Alignment Process
All the most recent alignments were produced by Richard Durbin's group at the Sanger
based on our analysis.sequence.index file which contains all our ILLUMINA sequence data
which has 70bp or longer reads. Older alignment releases also contain SOLID and 454 sequence
data and their alignment process is explained in the readme which sits with those
alignments
Illumina data was aligned with bwa v0.5.9 in 4 steps:
1. Index the reference fasta:
bwa index -a bwtsw $reference_fasta
2. For each fastq file:
bwa aln -q 15 -f $sai_file $reference_fasta $fastq_file
3. Create SAM files using bwa sampe or samse for paired-end or unpaired
reads respectively. For paired-end reads, the maximum insert size is
taken to be 3 times the expected insert size.
bwa sampe -a $max_insert_size -f $sam_file $reference_fasta $sai_files $fastq_files
bwa samse -f $sam_file $reference_fasta $sai_file $fastq_file
4. Create BAM from SAM, sort, fix mate information and add the MD tag:
samtools view -bSu $sam_file | samtools sort -n -o - samtools_nsort_tmp |
samtools fixmate /dev/stdin /dev/stdout | samtools sort -o - samtools_csort_tmp |
samtools fillmd -u - $reference_fasta > $fixed_bam_file
Bam Improvement
The run-level alignment BAMs are improved in various ways to help increase
the quality and speed of subsequent SNP calling that may be carried out on
them. For Illumina BAMs the following improvements were performed:
1. Reads undergo local realignment around known Indels using GATK
IndelRealigner.
java $jvm_args -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R $reference_fasta -o $intervals_file -known $known_indels_file(s)
java $jvm_args -jar GenomeAnalysisTK.jar -T IndelRealigner -R $reference_fasta -I $bam_file -o $realigned_bam_file -targetIntervals $intervals_file -known $known_indels_file(s) -LOD 0.4 -model KNOWNS_ONLY -compress 0 --disable_bam_indexing
where the known Indels are from the following:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.indels_mills_devine_hg19_leftAligned_collapsed_double_hit.indels.sites.vcf.gz
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.low_coverage_vqsr.20101123.indels.sites.vcf.gz
2. Realigned BAMs have their read qualities recalibrated with GATK
CountCovariates and TableRecalibration.
java $jvm_args -jar GenomeAnalysisTK.jar -T CountCovariates -R $reference_fasta -I $realigned_bam_file -recalFile recal_data.csv -knownSites $known_sites_file(s) -l INFO -L '1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;X;Y;MT' -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate
java $jvm_args -jar GenomeAnalysisTK.jar -T TableRecalibration -R $reference_fasta -recalFile recal_data.csv -I $realigned_bam_file -o $recalibrated_bam_file -l INFO -compress 0 --disable_bam_indexing
where the known sites for recalibration are from dbSNP135:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_mapping_resources/ALL.wgs.dbsnp.build135.snps.sites.vcf.gz
3. samtools calmd is applied to the recalibrated BAMs, which fixes the NM
tags and introduces BQ tags which can be used during SNP calling.
samtools calmd -Erb $recalibrated_bam_file $reference_fasta > $bq_bam_file
For Exome Solid BAMs, only step 2 was performed.
For low coverage Solid BAMs, only steps 1 and 2 were performed.
Release bam file production
The improved BAMs are merged together to get the release BAM files available
for download. Release BAM files therefore contain reads from multiple
readgroups. Production proceeds broadly as follows:
1. Run-level BAMs have extraneous tags (OQ, XM, XG, XO) stripped from them,
to reduce total file size by around 30%.
2. Tag-stripped BAMs are merged to the library level with Picard. (merging
to a given level means to create a single BAM file containing all the
readgroups that share a given member of that level; in this case it means
a BAM will be made for each library, each library BAM containing all the
reads from the run-level BAMs associated with that library).
java $jvm_args -jar MergeSamFiles.jar INPUT=$tag_stripped_bam_file(s) OUTPUT=$library_level_bam VALIDATION_STRINGENCY=SILENT
3. PCR duplicates are marked in library-level BAMs using Picard
MarkDuplicates.
java $jvm_args -jar MarkDuplicates.jar INPUT=$library_level_bam OUTPUT=$markdup_bam_file ASSUME_SORTED=TRUE METRICS_FILE=/dev/null VALIDATION_STRINGENCY=SILENT
4. Duplicate-marked library-level BAMs are merged to the platform level.
java $jvm_args -jar MergeSamFiles.jar INPUT=$markdup_bam_file(s) OUTPUT=$platform_level_bam VALIDATION_STRINGENCY=SILENT
5. Platform-level BAMs are split into chromosomes (and other, see above)
BAMs, which are the release files.
C. QA of the alignment data by DCC
Platform-level bams for individuals are checked by a QA process by DCC
before they are released to the ftp site. Here lists the QA measurements and criteria:
1. Check md5 for each file to make sure the file was not corrupted during transfer.
2. Check coverage status. For exome runs this is done by runnung calculate_HsMetrics
(a function in PICARD) to evaluate EXOME sequence coverage. Samples with less than
70% 20x coverage in targetted regions are considered coverage too low. The bait file
used in this calculation is ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/exome_pull_down_targets_phases1_and_2/20120518.consensus.annotation.bed
For low coverage samples we use the bas files supplied with the bams. These contain
the number of mapped bases and the number of duplicate bases. We calculate the non
duplicated mapped coverage. This must be 3x or greater (consdering a genome size of 2.75Gb)
3. Check for excessive number of short insertions or deletions (<=6bp). If the ratio of
short insertion counts to short deletion counts is greater than 5, the sample fails.
The reverse ratio of short deletion counts to short insertion counts is also checked
the same way. (This uses a custom script, if you would like a copy please email info@1000genomes.org)
4. Check for sample contamination using VerifyBAMId developed by Hyun Kang in University of Michigan.
VerifyBamId is a software that verifies whether the reads in particular file match previously
known genotypes for an individual (or group of individuals), and checks whether the reads are
contaminated as a mixture of two samples. verifyBamID can detect sample contamination and swaps
when external genotypes are available. When external genotypes are not available, verifyBamID still
robustly detects sample swaps. Genotype data for most phase1 and 2 samples typed on OMNI chips and
genotype data for most phase3 samples typed on Affymetrix chips were used in this QA process. The original OMNI
genotype vcf files can be found at /nfs/1000g-archive/vol1/ftp/technical/working/20120131_omni_genotypes_and_intensities/Omni25_genotypes_2141_samples.b37.vcf.gz
and Affy genotype data can be found at /nfs/1000g-archive/vol1/ftp/technical/working/20121128_corriel_p3_sample_genotypes/. We divided the vcf files by sample
population when ran VerifyBamId.
The CHIP_MIX and FREE_MIX cutoff ranges from 2-3.5% for different BAM releases.
http://genome.sph.umich.edu/wiki/VerifyBamID
5. Check for completeness of data against sequence index:
a. All reads and runs for an individual should be included in
corresponding BAM files
b. All reads and runs from a BAM file should be from the same
individual
It is important to note that the 20130502 bam release was unusual as rather than all QC passed bams
did not make it into the data/XXXXXX/alignment directories
For the final round of analysis for the project we only used ILLUMINA sequence data which was 70bp or greater.
The group decided for the final variant calling only samples with both low coverage and exome sequence would
be considered. This means there were a small number of samples with just exome or just low coverage alignments.
The sequence data is still included in the sequence.index file and while the alignments are not being used in
our variant calling they are still available to be downloaded.
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/phase3_EX_or_LC_only_alignment/
These bam files have their own alignment index which you can find in this directory.
D. bas files
.bas files are bam statistic files, one line per readgroup, columns separated by
tabs. The first line is a header that describes each column. The first 6 columns
provide meta information about each readgroup.
The remaining columns provide various statistics about the readgroup, calculated
by going through the release bams. Where data isn't available to calculate the
result for a column the default value will be 0.
Each column is described in detail below.
Column 1 'bam_filename': The DCC bam file name in which the readgroup data can
be found.
Column 2 'md5': The md5 checksum of the bam file named in column 1.
Column 3 'study': The SRA study id this readgroup belongs to.
Column 4 'sample': The sample (individual) identifier the readgroup came from.
Column 5 'platform': The sequencing platform (technology) used to sequence the
readgroup.
Column 6 'library': The name of the library used for the readgroup.
Column 7 'readgroup': The readgroup identifier. This is unique per .bas file.
The remaining columns summarise data for reads with this RG tag in the bam
file given in column 1.
Column 8 '#_total_bases': The sum of the length of all reads in this readgroup.
Column 9 '#_mapped_bases': The sum of the length of all reads in this readgroup
that did not have flag 4 (== unmapped).
Column 10 '#_total_reads': The total number of reads in this readgroup.
Column 11 '#_mapped_reads': The total number of reads in this readgroup that did
not have flag 4 (== unmapped).
Column 12 '#_mapped_reads_paired_in_sequencing': As for column 10, but also
requiring flag 1 (== reads paired in sequecing).
Column 13 '#_mapped_reads_properly_paired': As for column 10, but also requiring
flag 2 (== mapped in a proper pair, inferred during alignment).
Column 14 '%_of_mismatched_bases': Calculated by summing the read lengths of all
reads in this readgroup that have an NM tag, summing the edit distances
obtained from the NM tags, and getting the percentage of the latter out of
the former to 2 decimal places.
Column 15 'average_quality_of_mapped_bases': The mean of all the base qualities
of the bases counted for column 8, to 2 decimal places.
Column 16 'mean_insert_size': The mean of all insert sizes (ISIZE field) greater
than 0 for properly paired reads (as counted in column 12) and with a
mapping quality (MAPQ field) greater than 0. Rounded to the nearest whole
number.
Column 17 'insert_size_sd': The standard deviation from the mean of insert sizes
considered for column 15. To 2 decimal places.
Column 18 'median_insert_size': The median insert size, using the same set of
insert sizes considered for column 15.
Column 19 'insert_size_median_absolute_deviation': The median absolute deviation
of the column 17 data.
Column 20 '#_duplicate_reads': The number of reads which were marked as
duplicates
Column 21' #_duplicate_bases': The number of bases which were narked as
duplicated