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

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

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

Подготовительные операции.

Для начала объединим полученные в предыдущей части ped-файлы в  общий набор данных. Задача выполняется тривиальной командой —bmerge —make-recode.

Затем добавляем в набор данных подготовленные ранее файлы, содержащие в себе данные о снип-вариантах, обнаруженных в аутосомных хромосомах «тирольского человек Этци». Мимоходом отмечу, что оные данные были извлечены из генома Этци благодаря использованию несколько иной программы (GATK) и методологии, поэтому я не использовал их в предыдущей части туториала.

 

----------------------------------------------------------@
| PLINK! | v1.07 | 10/Aug/2009 |
|----------------------------------------------------------|
| (C) 2009 Shaun Purcell, GNU General Public License, v2 |
|----------------------------------------------------------|
| For documentation, citation & bug-report instructions: |
| http://pngu.mgh.harvard.edu/purcell/plink/ |
@----------------------------------------------------------@
Skipping web check... [ --noweb ] 
Writing this text to log file [ lastreference.log ]
Analysis started: Wed Apr 10 16:37:28 2013
Options in effect:
 --noweb
 --bfile lastreference
 --bmerge Otzi2.bed Otzi2.bim Otzi2.fam
 --make-bed
 --out lastreference
Reading map (extended format) from [ lastreference.bim ] 
424693 markers to be included from [ lastreference.bim ]
Reading pedigree information from [ lastreference.fam ] 
2618 individuals read from [ lastreference.fam ] 
0 individuals with nonmissing phenotypes
Assuming a disease phenotype (1=unaff, 2=aff, 0=miss)
Missing phenotype value is also -9
0 cases, 0 controls and 2618 missing
1388 males, 759 females, and 471 of unspecified sex
Warning, found 471 individuals with ambiguous sex codes
Writing list of these individuals to [ lastreference.nosex ]
Reading genotype bitfile from [ lastreference.bed ] 
Detected that binary PED file is v1.00 SNP-major mode
Using merge mode 1 : consensus call (default)
138056 markers to be merged from [ Otzi2.bim ]
Of these, 0 are new, 138056 already exist in current data
1 individuals merged from [ Otzi2.fam ] 
Of these, 1 were new, 0 were already in current data
Detected that binary PED file is v1.00 SNP-major mode
0 individuals with nonmissing phenotypes
Assuming a disease phenotype (1=unaff, 2=aff, 0=miss)
Missing phenotype value is also -9
0 cases and 0 controls
Before frequency and genotyping pruning, there are 424693 SNPs
2619 founders and 0 non-founders found
197 SNPs with no founder genotypes observed
Warning, MAF set to 0 for these SNPs (see --nonfounders)
Writing list of these SNPs to [ lastreference.nof ]
Total genotyping rate in remaining individuals is 0.587873
0 SNPs failed missingness test ( GENO > 1 )
0 SNPs failed frequency test ( MAF < 0 )
After frequency and genotyping pruning, there are 424693 SNPs
After filtering, 0 cases, 0 controls and 2619 missing
After filtering, 1388 males, 759 females, and 472 of unspecified sex
Writing pedigree information to [ lastreference.fam ] 
Writing map (extended format) information to [ lastreference.bim ] 
Writing genotype bitfile to [ lastreference.bed ] 
Using (default) SNP-major mode
Analysis finished: Wed Apr 10 16:40:14 2013

Решение проблемы совместимости данных, хранящихся в разных геномных билдах (сборках).

Одной из cамых болезненных проблем, c которыми сталкивается исследователь геномного разнообразия, является проблема совместимости геномных билдов (или сборок).  Определенные затруднения вызваны тем, что данные по разным популяциям могут быть приведены в разных билдах. Более ранние выборки были представлены в NCBI билде hg/b36, однако в последнее время все больше новых данных поставляется в новом NCBI билде b37. Во-многом это объясняется успехом проекта 1000 геномов, в котором как раз и используются мэппинг (генетических и физических позиций) NCBI билда b37. Наверное, по этой же причине коммерческие компании типа 23andme и FTDNA также сменили версию билда на b37.

Мы также столкнулись с аналогичной проблемой, поскольку снип-данные были извлечены из древней ДНК  с использованием более старого билда.  При решении этой проблемы мы руководились советами авторов Enigma Cookbook 1000 Genomes Imputation.

В практическом плане имеется целый ряд различий между двумя сборками. Были выявлены три основных типа проблем в build36, которые решены в build37:

● были найдены  SNP-ы встречающиеся два раза под разными rs-индексами
● обнаружено несколько SNP-ов, отображенных на неверной хромосому
● были найдены многочисленные SNP-ы c неверными физическими координатами (но на верной
хромосомы)
Для преобразования данных из сборки build36 в сборку build37  нам необходимо переназначить SNP-позиции
позиций, перечисленных в 1KGP, используя следующий код.  Существует небольшое количество полиморфизмов (около 312) в референсной панели 1000 геномов (1KGP), которые отображаются на БОЛЕЕ чем одной хромосоме, поэтому перед тем как продолжить, нам необходимо отбросить эти снипы.

Референсная панель (data for 41 million markers, of which ~23 million are monomorphic in Caucasians) скачивается следущей последовательностью команд:

mkdir 1KGPref
cd 1KGPref
wget "http://enigma.loni.ucla.edu/wp-content/uploads/2012/07/
v3.20101123.ENIGMA2.EUR.20120719.vcf.tgz"
wget "http://enigma.loni.ucla.edu/wp-content/uploads/2012/07/
v3.20101123.ENIGMA2.EUR.20120719.extras.tgz"
tar -zxvf v3.20101123.ENIGMA2.EUR.20120719.vcf.tgz
tar -zxvf v3.20101123.ENIGMA2.EUR.20120719.extras.tgz

 

Следующий код загружает генетическую карту панели 1KGP и создает список, который будет использоваться для фильтрации генотипированных данных перед фазированием:

## For the following commands in green use the clean “lastQC2” files
if you had to remove duplicate markers
#Join the genotyped bim file with the reference allele lists
##reformat the lastQC.bim file
awk '{print $2,$1,$3,$4,$5,$6}' lastQC.bim > tempQC.bim
##Join the two files
awk 'NR==FNR{s=$1;a[s]=$0;next} a[$1]{print $0 " "a[$1]}' tempQC.bim
1kgp.alleles > merged.alleles
## selects SNPS showing different alleles in the two files
awk '{ if ($2!=$8 && $2!=$9) print $1}' merged.alleles > flip.list
plink --bfile lastQC --extract 1kgp.snps --update-map 1kgp.chr --
update-chr --flip flip.list --make-bed --out temp --noweb
plink --bfile temp --update-map 1kgp.bp --geno 0.05 --mind 0.05 --
make-bed --out lastQCb37 --noweb
wc -l lastQCb37.bim
## Make list of males and females for writing out the X chromosome
awk '{ if($5==1) print $1, $2}' lastQCb37.fam > male.list
awk '{ if($5==2) print $1, $2}' lastQCb37.fam > female.list
##Check that your dataset is properly split by gender by opening
male.list and female.list in a text editor. Also check that total
numbers make sense.
wc -l female.list
wc -l male.list
@----------------------------------------------------------@
| PLINK! | v1.07 | 10/Aug/2009 |
|----------------------------------------------------------|
| (C) 2009 Shaun Purcell, GNU General Public License, v2 |
|----------------------------------------------------------|
| For documentation, citation & bug-report instructions: |
| http://pngu.mgh.harvard.edu/purcell/plink/ |
@----------------------------------------------------------@
Skipping web check... [ --noweb ] 
Writing this text to log file [ lastQCb37.log ]
Analysis started: Wed Apr 3 02:09:06 2013
Options in effect:
 --bfile temp
 --update-map 1kgp.bp
 --out lastQCb37
 --noweb
 --make-bed
Reading map (extended format) from [ temp.bim ] 
1032625 markers to be included from [ temp.bim ]
Reading pedigree information from [ temp.fam ] 
3606 individuals read from [ temp.fam ] 
4 individuals with nonmissing phenotypes
Assuming a disease phenotype (1=unaff, 2=aff, 0=miss)
Missing phenotype value is also -9
0 cases, 4 controls and 3602 missing
1741 males, 975 females, and 890 of unspecified sex
Warning, found 890 individuals with ambiguous sex codes
Writing list of these individuals to [ lastQCb37.nosex ]
Reading genotype bitfile from [ temp.bed ] 
Detected that binary PED file is v1.00 SNP-major mode
Reading new physical positions from [ 1kgp.bp ]
1032625 SNP positions read and updated
0 in data but not in [ 1kgp.bp ]
11797154 in [ 1kgp.bp ] but not in data
*** Implicit order changed from re-mapping ***
Before frequency and genotyping pruning, there are 1032625 SNPs
3606 founders and 0 non-founders found
273 heterozygous haploid genotypes; set to missing
Writing list of heterozygous haploid genotypes to [ lastQCb37.hh ]
Total genotyping rate in remaining individuals is 0.309747
0 SNPs failed missingness test ( GENO > 1 )
0 SNPs failed frequency test ( MAF < 0 )
After frequency and genotyping pruning, there are 1032625 SNPs
After filtering, 0 cases, 4 controls and 3602 missing
After filtering, 1741 males, 975 females, and 890 of unspecified sex
Writing pedigree information to [ lastQCb37.fam ] 
Writing map (extended format) information to [ lastQCb37.bim ] 
Writing genotype bitfile to [ lastQCb37.bed ] 
Using (default) SNP-major mode
Analysis finished: Wed Apr 3 02:13:38 2013

Теперь, когда координаты и странды ДНК приведены в порядок, мы можем соединять данные древней ДНК с нашей главной референсной панелью.

В следующей части мы покажем конкретные результаты сравнения древней и современной ДНК.