4.6 liftOver
Lifting GERP scores from one reference to another
liftOver allows for a simple way to lift over conservation GERP scores from one reference to another. Multiple ATLAS tasks, such as mutationLoad and estimateErrors, recommend using information at conserved sites. The liftOver task offers an easy solution to identify these sites in reference genomes, for which GERP scores are not already available.
This task requires four steps (as shown in the example below).
- Downloading both the GERP scores and the matching reference genome of one of the available species. Choose the species most closely related to the target reference coordinate system of choice.
Bed2Fastqmode: Using the downloaded GERP scores and the matching reference genome to create fastq reads. The task will create one fastq read for every GERP position, including flanking bases on both sides of the bed position. The GERP score is encoded in the read name- Mapping: Mapping those fastq reads to the reference of interest
Bam2Bedmode: Creating the bed file of lifted over positions and their GERP scores in the reference coordinate system of choice
Note: this method has proven to be accurate and productive if references are reasonably close (eg. tested with goat GERP scores lifted over to ibex or with horse GERP scores lifted over to hippo or with one Darwin’s finch GERP scores lifted over to another) and works best with conserved sites. Seek out the available GERP scores of a species as closely related as possible to the target species and constrain the liftOver to positions with conserved GERP scores.
Extensive tutorial coming soon!
4.6.1 Bed2Fastq mode
4.6.1.1 Bed2Fastq input
Required inputs:
--fasta original.fasta |
Input reference genome, for which GERP scores are available |
--bed gerp.bed |
Input bed file with GERP scores in the reference coordinate system of original.fasta |
Example bed file containing GERP scores:
chr1 1000 1001 7.2
chr1 1001 1002 9
chr2 3000 3001 5Specific Parameters:
--flank length |
to provide length of the flanking bases on both sides of the bed position (Default = 100) |
4.6.2 Bam2Bed mode
4.6.3 Example of the workflow
#! /bin/bash
# Set atlas path
atlas=$(dirname "$0")/../build/atlas
## PREPARATION
# Simulate a reference genome
$atlas simulate --chrLength 100000
# Create an example bed file of how downloaded GERP scores might look
# To simplify for this example, we create a bed file with every 200th position and a random GERP score. Be aware that this is not a realistic example.
printf "chr1\t0\t1\t7\n" > GERP.bed
for i in {200..100000..200}; do
printf "chr1\t$i\t$((i+1))\t$((RANDOM % 10))\n" >> GERP.bed
done
# Create a related reference for which GERP scores should be lifted over from the original reference. In a realistic scenario, this reference would be the different species of interest, for which no GERP scores are available.
# Be aware that we adjust the simulated reference only as a illustration of how to use the liftOver task with two related references.
sed '50d' ATLAS_simulations.fasta > target_reference.fasta
## STEP1: Create fastq reads at all positions in the given GERP.bed file using the original reference and the liftOver Bed2Fastq mode.
$atlas liftOver --mode Bed2Fastq --bed GERP.bed --fasta ATLAS_simulations.fasta
## STEP2: Map the fastq reads to the related reference of interest using your mapper of choice
bwa index target_reference.fasta # indexing the reference
bwa mem -M -t 1 target_reference.fasta ATLAS_Bed2Fastq.fastq | samtools view -bSh - | samtools sort -o ATLAS_liftOver.bam -
## STEP3: Create the final bed file of lifted over positions including their GERP scores by using the bam file from step 2 in the liftOver Bam2Bed mode
$atlas liftOver --mode Bam2Bed --bam ATLAS_liftOver.bam