Практические рекомендации по работе с данными древней ДНК

Опубликовано автором verenich

В отличие от большинства записей в моем блоге,  эта запись будет посвящена практическим аспектам работы с геномными данными доисторических останков. В целях экономии времени и пространства, я пока не буду затрагивать вопросы связанные с чисто технической стороной работы с древней ДНК, тем более что ответы на эти вопросы неплохо освящены в соответствующей литературы(кратких конспектов).

Следует также заметить, что стиль изложения материала в данной заметке намерено упрощен в целях облегчения материала. Исходя из этого следует помнить, что чтение этого материала никоим образом не заменит собой более тщательного и глубокого ознакомления с исследовательской методологией.

В качестве примера в нашем туториале мы будем использовать данные, любезно предоставленные авторами работы P Skoglund, H Malmström, M Raghavan, J Storå, P Hall, E Willerslev, MTP Gilbert, A Götherström* & M Jakobsson* (2012) Origins and genetic legacy of Neolithic farmers and hunter-gatherers in Europe, и данные работы  Federico Sánchez-Quinto, Hannes Schroeder, Oscar Ramirez, María C. Ávila-Arcos, Marc Pybus, Iñigo Olalde, Amhed M.V. Velazquez, María Encina Prada Marcos, Julio Manuel Vidal Encinas, Jaume Bertranpetit, Ludovic Orlando, M. Thomas P. Gilbert, Carles Lalueza-Fox Genomic Affinities of Two 7,000-Year-Old Iberian Hunter-Gatherers.

Для успешного прохождения туториала нам потребуется:
1) наличие мотивации и желание изучить основы практической геномики

2) посколько большинство инструментов задействованных в данном туториале написны под Unix, то необходимо наличие опыта работы с Unix shell: желательно также иметь доступ к значительным вычислительным мощностями (некоторые из операций описанных ниже я производил в вычислительном кластере Тартуского университета).

3) пакет samtools последней версии

4) пакет snpEFF/snpSift

5) пакет vcftools и программа Plink

6) FASTA-файл с человеческий референсным геномом в  версии билда hg18: я рекомендую использовать модифицированную версию файла, в котором старый референсный митосиквенс заменен на новый референсный митосиквенс (rCRS:NC_012920 gi:251831106).

7) каталог генетических полиморфизмов dbSNP в версии билда hg18.

I этап — bamtools.

Используемые в нашем туториале исходные файлы представлены в формате bam — бинарном варианте стандартного файла SAM используемого для хранения элайнментов сиквенсов.
В нашем случае исходные файлы представляют собой конечный продукт,  в котором уже удалены дупликаты и артифакты клонирования в ходе PCR. Поэтому мы можем сразу же приступить к следущему этапу — объединения файлов bam в один общий файл:

samtools merge AjvIre.bam SNPs_Ajv52_r1_hits_rmdup.bam SNPs_Ajv52_r2_hits_rmdup.bam SNPs_Ajv70_r1_hits_rmdup.bam SNPs_Ajv70_r2_hits_rmdup.bam SNPs_Ire8_r1_hits_rmdup.bam SNPs_Ire8_r2_hits_rmdup.bam  
samtools merge GokSte.bam SNPs_Ste7_r1_hits_rmdup.bam SNPs_Ste7_r2_hits_rmdup.bam SNPs_Ste7_r3_hits_rmdup.bam SNPs_Ste7_r4_hits_rmdup.bam\

Далее, мы провидем сортировку файлов по контигам (контиг — это набор упорядоченных перекрывающихся клонов ДНК, охватывающих всю хромосому ил и какой-либо ее участок):

 samtools sort AjvIre.bam AjvIre.sorted.bam
 samtools sort GokSte.bam GokSte.sorted.bam
 samtools sort BRA.bam BRA.sorted.bam

Скачиваем референсный файл билда hg18

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg18/bigZips/hg18.2bit
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
./twoBitToFa hg18.2bit hg18.fa

Производим индексацию референсного файла человеческого генома (билд hg18) и сравниваем систему обозначения хромосомных контигов с аналогичной системой в наших образцах древних геномов:

samtools faidx hg18.fa
chr10    135374737    7    50    51
chr10_random    113275    138082253    50    51
chr11    134452384    138197801    50    51
chr11_random    215294    275339247    50    51
chr12    132349534    275558854    50    51
chr13    114142980    410555386    50    51
chr13_random    186858    526981240    50    51
chr14    106368585    527171843    50    51
chr15    100338915    635667807    50    51
chr15_random    784346    738013515    50    51
chr16    88827254    738813555    50    51
chr16_random    105485    829417369    50    51
chr17    78774742    829524971    50    51
chr17_random    2617613    909875222    50    51
chr18    76117153    912545195    50    51
chr18_random    4262    990184706    50    51
chr19    63811651    990189061    50    51
chr19_random    301858    1055276960    50    51
chr1    247249719    1055584862    50    51
chr1_random    1663265    1307779589    50    51
chr20    62435964    1309476127    50    51
chr21    46944323    1373160818    50    51
chr21_random    1679693    1421044042    50    51
chr22    49691432    1422757336    50    51
chr22_random    257318    1473442611    50    51
chr22_h2_hap1    63661    1473705091    50    51
chr2    242951149    1473770032    50    51
chr2_random    185571    1721580217    50    51
chr3    199501827    1721769506    50    51
chr3_random    749256    1925261383    50    51
chr4    191273063    1926025631    50    51
chr4_random    842648    2121124169    50    51
chr5    180857866    2121983676    50    51
chr5_random    143687    2306458713    50    51
chr5_h2_hap1    1794870    2306605288    50    51
chr6    170899992    2308436062    50    51
chr6_random    1875562    2482754067    50    51
chr6_cox_hap1    4731698    2484667156    50    51
chr6_qbl_hap2    4565931    2489493503    50    51
chr7    158821424    2494150759    50    51
chr7_random    549659    2656148625    50    51
chr8    146274826    2656709284    50    51
chr8_random    943810    2805909620    50    51
chr9    140273252    2806872313    50    51
chr9_random    1146434    2949951044    50    51
chrM    16571    2951120413    50    51
chrX    154913754    2951137322    50    51
chrX_random    1719168    3109149365    50    51
chrY    57772954    3110902923    50    51

samtools view -H AjvIre.sorted.bam

@HD    VN:1.0    SO:unsorted@PG    ID:dvtgmlqtca    PN:stampy    VN:1.0.10_(r854)    CL:-g hg18 -h hg18 --solexa --sensitive -f sam -o output/stampy_Ajv52_r1_aln1.sam -M /bubo/proj/b2010050/private/seqdata/neolitikum/Neolitisar/pruned/Ajv52_r1_trimmed.txt@CO    TM:Tue, 30 Nov 2010 12:35:43 CET    WD:/bubo/proj/b2010050/private/program/stampy-1.0.10    HN:q207.uppmax.uu.se    UN:pontuss@SQ    SN:NC_000001.9    LN:247249719    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000002.10    LN:242951149    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000003.10    LN:199501827    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000004.10    LN:191273063    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000005.8    LN:180857866    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000006.10    LN:170899992    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000007.12    LN:158821424    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000008.9    LN:146274826    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000009.10    LN:140273252    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000010.9    LN:135374737    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000011.8    LN:134452384    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000012.10    LN:132349534    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000013.9    LN:114142980    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000014.7    LN:106368585    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000015.8    LN:100338915    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000016.8    LN:88827254    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000017.9    LN:78774742    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000018.8    LN:76117153    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000019.8    LN:63811651    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000020.9    LN:62435964    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000021.7    LN:46944323    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000022.9    LN:49691432    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000023.9    LN:154913754    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_000024.8    LN:57772954    AS:hg18_ncbi36_rCRS    SP:human@SQ    SN:NC_012920.1    LN:16560    AS:hg18_ncbi36_rCRS    SP:human

Итак,  при сравнении вышеупомянутых двух файлов, мы видим что обозначение в референсном генома  отличается от обозначения начала хромосом в файле AjvIre (вместо традиционного обозначения chr1…chrM в этом файле используется номер сиквенса хромосомы в Генбанке, например сиквенс первой хромосомы — NC_000001.9, и т.д.).
Эта проблема решается сравнительно легко с помощью редактирования заголовка bam файла (заменой SO:unsorted@PG  на SO:sorted@PG и номеров Генбанка на порядковый номер хромосом) и следущих комбинаций директив samtools:

samtools view -H AjvIre.sorted.bam > originalheader
gedit originalheader
samtools reheader newheader AjvIre.reh.sorted.bam

Аналогичные операции производим и с файлом GokSte.sorted.bam. Файл Bra.sorted.bam редактировать нет надобности, поскольку обозначения хромосом соответствуют обозначению хромосом в референсном файле.

Таким образом, после выполнения означенных выше операций, мы подошли к самой важной процедуре — snp and indel calling, то есть «вызову» (определению) снипов и инделов в наших отсортированных и модифицированных bam файлах.

Нужно сразу отметить, что процедура нахождения генетических вариантов в древней ДНК существенно отличается от аналогичной процедуры в случае с современной ДНК.  Поэтому приходится применять фильтры samtools, которые в большинстве рутинных анализов просто не используются. Я не буду объяснять, что означает каждый из используемых фильтров. Достаточно будет сказать, что я следую рекомендациям профессора Понтуса Скоглунда.  Принимая во внимание ресурсоемкость операции нахождения генетических вариантов, я задействовал возможности тартуского вычислительного центра (ниже приведен пример с BRA.srt.bam):

qsub runSamstools.sh

#!/bin/bash
# This file is runSamtools
#
#PBS -N Samtools
#PBS -m be
#PBS -k oe
#PBS -l walltime=01:30:00
#PBS -l nodes=4:ppn=8
#PBS -l vmem=4gb
#PBS -d /storage/hpchome/vadim78

cd /storage/hpchome/vadim78/conversion/ancient
module load storage_software
samtools mpileup BRA.srt.bam -q 30 -Q 15 -uf hg18.fa |
/storage/hpchome/vadim78/samtools/bcftools/bcftools view -vcg - > BRA.vcf

II. Аннотация VCF файлов — snpSift.

Итак,  мы получили три файла VCF, которые содержат в себе информацию о найденных генетических вариантах — инделах и снипах.  При визуальном осмотре файлов сразу же бросается в глаза отсутствие идентификаторов снипов/инделов. Вместо привычных rs-id, варианты индексированы с помощью точек . Поскольку нам необходима для дальнейшего анализа традиционная система обозначения, мы должны произвести аннотирование файлов. Путем метода проб и ошибок я выбрал самую удобный для начинающих геномиков пакет snpEff.

В качестве источника аннотирования мы используем каталог генетических вариантов dbSNP, который содержит не только rs-индексы снипов, но и широкий спектр данных о функциональных связах снипа, в том числе и о генетических ассоциациях. Но мы не будем рассматривать весь спектр данных, поскольку нас интересуют индексы снипов.

Аннотирование  индексов снипов  в VCF-файлах выполняется с помощью  несложной командой (ниже приведен пример командной строки для файла GokSte.vcf).

java -Xmx2g -jar SnpSift.jar annotateMEM -id  ../dbsnp.vcf ../GokSte.vcf > GokSte.annotate.vcf

Очевидно, что древняя ДНК содержит значительное число новых истинных и ложных снипов, которых нет в индексах dbSNP.  В нашем туториале мы ограничимся лишь известными снипами, и поэтому отфильтруем «новельные» снипы.

java -Xmx2g -jar SnpSift.jar filter "(exists ID) & ( ID =~ 'rs' )"  GokSte.annotate.vcf > GokSte.snp.vcf

III.  Фильтрация снипов в vcftools

Как я указывал в предыдущем разделе, файлы VCF содержат в себе информацию о всех найденных генетических вариантах — инделах и снипах.  Несмотря на всю важность инделов в определение вариативности генофонда популяций,  во многих популяционно-генетических исследованиях явное предпочтение отдается снипам. Принимая это во внимание, я решил отсеять инделы в сторону и трансформировать файл VCF в более традиционный формат Plink PED:

./vcftools/bin/vcftools --vcf  GokSte.snp.vcf --remove-indels  --plink --out GokSte

На выходе, мы получил файл Plink PED, о котором мы поговорим в следущей части туториала.

Продолжение следует.