Как объединить данные SNP с эталонным геномом?

Мои данные

У меня есть файл 23andMe со списком SNP в виде:

rsid chromosome position genotype rsXXXXX 1 PPPPPP CT rsXXXXX 1 PPPPPP GG

Поля разделены TAB, и каждая строка соответствует одному SNP. Для каждого SNP предоставляются четыре поля данных.

  1. Идентификатор (rsid или внутренний идентификатор)
  2. Его расположение на эталонном геноме.
    • Хромосома, на которой он расположен.
    • Положение внутри хромосомы находится на .
  3. Вызов генотипа ориентирован относительно плюсовой цепи референсной последовательности человека.

Эталонный геном представляет собой сборку 37 сборки человека (также известную как выпуск аннотаций 104).

Мой вопрос

Как мне объединить SNP в эталонный геном?

Например, возьмем первую строку в моем файле SNP:

rsXXXXX 1 PPPPPP CT

Часть 1

Я вижу, что мне нужно заменить нуклеотид в положении PPPPPP на хромосоме 1 эталонного генома на нуклеотид из поля генотипа, но какой нуклеотид я должен использовать? С или Т? И почему?

Часть 2

С чего я должен начать отсчет эталонного генома? Глядя на хромосому 1 сборки человеческой сборки 37, первые ~ 10 000 символов (исключая описание в первой строке) — это N. Является ли первый N номером 1? например. Если бы PPPPPP было равно 100 000, заменил бы я 100 000-й символ в эталонном геноме правильным нуклеотидом из части 1 этого вопроса? Или я должен начать считать с первого символа, отличного от N, в файле fasta?

Ответы (3)

Во-первых, вам нужно знать, к какой последовательности генома относится файл SNP. Они, должно быть, упомянули референсную последовательность, которую они использовали.

Как уже упоминалось, случай CTгетерозиготности. Если вы просто хотите отметить изменения, отбросьте остаток, который уже присутствует в эталонном геноме, и используйте другой аллель. Тем не менее, вы хотите отслеживать гаплотип, тогда вам нужно убедиться, что набор SNP происходит из одной и той же хроматиды. Это сложно - вы все еще можете узнать это для SNP, которые достаточно близки, чтобы их можно было отобразить одним считыванием, но это почти невозможно для SNP, которые достаточно хорошо разделены.

Как сказал Эндре, нужно начинать с первого нуклеотида. Тем не менее, кажется сомнительным, что вы получаете   ( Н Н Н Н ) н в начале хромосомы 1. Полно собранные хромосомы таких участков не имеют. Ниже приведены первые 10 строк файла fasta хромосомы 1. Посмотреть на себя.

>gi|568815364|ref|NT_077402.3| Homo sapiens chromosome 1 genomic scaffold, GRCh38 Primary Assembly HSCHR1_CTG1
TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC
CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAACCCTAACCCTAA
CCCTAACCCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACCCTAACCCTAACCCTAACCCT
AACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAAACCCTAAACCCTAACCCTAACCCTAACCCTA
ACCCTAACCCCAACCCCAACCCCAACCCCAACCCCAACCCCAACCCTAACCCCTAACCCTAACCCTAACC
CTACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCCTAACCCTAACCCTAACCCTA
ACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACCCTCGCGGTACCCTCAGCCGGCCCG
CCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTACCACCGAAATCTGTGCAGAGGAC
AACGCAGCTCCGCCCTCGCGGTGCTCTCCGGGTCTGTGCTGAGGAGAACGCAACTCCGCCGTTGCAAAGG

Как заменить Н т час Остаток - довольно простая задача. Но это вопрос программирования, а не предмет обсуждения на этом форуме. Предполагая, что вы решили проблему части 1 и имеете отсортированный файл, разделенный табуляцией, например:

chromosome  position    residue
 1           79989           G
 1           100232          T
 3           341342          A

Этот скрипт может быть не самым лучшим, но он будет работать в терминале linux/*nix/Cygwin, чтобы заменить остатки (убедитесь, что у вас есть gawk version >=4.0):

gawk -F "\t" '(FNR==1){x++} (x==1){a[$1][$2]=$3;next} (x==2){if($0~/>/){h=$0;sub(/^.*chromosome /,"",h);sub(/ .*/,"",h)} else{seq[h]=seq[h]$0}} END{for(i in a){s=0; for(j in a[i]){m=m substr(seq[i],s,j-1) a[i][j];s=j+1} m=m substr(seq[i],s); print ">Chr"i"\n"m}}' SNP_file Genome.fa | fold -w 60
К сожалению, я не могу объяснить, как работает скрипт, на этом форуме
"TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC" до позиции 10 000 согласно UCSC. И посмотрите на это, как я уже сказал, это повторяющееся и ужасное. OP просто имеет замаскированную последовательность.
Хм.. да.. не заметил.. тем не менее начало с остатка №1..

Генетика 101, у вас есть 2 копии всей вашей ДНК в каждой позиции, одна копия от вашей матери, одна от вашего отца. Итак, для «CT» у вас есть одна копия с буквой C и одна с буквой T.

И да, первые несколько тысяч или миллионов букв — это нормально. Геном здесь повторяющийся и неприятный, но он все равно учитывается для целей нумерации.

Честно говоря, я бы не стал делать это с гигантским текстовым файлом генома. Просто найдите свой SNP на ensembl.org, используя номер rs, и вы получите SNP, некоторую фланкирующую последовательность и некоторый контекст. Поищите его в PubMed, если хотите узнать, появлялся ли он когда-либо в какой-либо публикации.

Часть 1:

По словам Лиора Пахтера , данные 23andme не поэтапны. Это означает, что для каждой записи в поле генотипа вы не знаете, из какой копии хромосомы она произошла. Это происходит потому, что современные платформы микрочипов не могут определить, из какой из двух копий хромосомы произошел snp.

Вы можете решить эту проблему для большинства snps, сравнив ваши аллели с эталонным геномом, но это потребует некоторых навыков программирования. Вы можете использовать https://github.com/endrebak/qc_gwas в качестве примера, который делает то же самое, но для файлов plink.

Часть 2:

Я предполагаю, что вы хотели бы сделать это программно, а не путем копирования и вставки snps в эталонный геном.

Короткий ответ: первый N — это первый нуклеотид. Но вам лучше использовать пакет, такой как Biopython , чтобы сделать подсчет за вас, он может быть более запутанным, чем вы думаете (например, вам нужно настроить окончания строк в файле fasta).

Как объединить данные SNP с эталонным геномом?
СпросиСетьПолитика конфиденциальности1y, 0v