Этот пример показывает, как сравнивать целые геномы для организмов, что позволяет сравнивать организмы с очень разным разрешением относительно сравнений отдельных генов. Вместо того, чтобы просто сосредоточиться на различиях между гомологичными генами, вы можете получить представление о крупномасштабных особенностях геномной эволюции.
В этом примере используются два штамма Chlamydia, Chlamydia trachomatis и Chlamydophila pneumoniae. Это близкородственные бактерии, которые вызывают разные, хотя и очень распространенные, заболевания у людей. Целые геномы доступны в базе данных GenBank ® для обоих организмов.
Вы можете загрузить эти геномы с помощью getgenbank функция. Сначала загрузите геном Chlamydia trachomatis. Обратите внимание, что геном круговой и длиной чуть более миллиона п.н. Эти последовательности достаточно велики, поэтому загрузка может занять некоторое время.
seqtrachomatis = getgenbank('NC_000117');
Далее загрузите Chlamydophila pneumoniae. Этот геном также круговой и немного длиннее при 1,2 Мбит/с.
seqpneumoniae = getgenbank('NC_002179');
Для удобства ранее загруженные последовательности включаются в MAT-файл. Следует отметить, что данные в публичных репозиториях часто обрабатываются и обновляются. Следовательно, результаты этого примера могут несколько отличаться при использовании актуальных наборов данных.
load('chlamydia.mat','seqtrachomatis','seqpneumoniae')
Очень простой подход для сравнения двух геномов заключается в выполнении парного выравнивания между всеми генами в геномах. Учитывая, что это бактериальные геномы, простым подходом было бы сравнение всех ORF в двух геномах. Однако данные GenBank включают больше информации о генах в последовательностях. Он хранится в поле CDS структуры данных. Chlamydia trachomatis имеет 895 кодирующих областей, в то время как Chlamydophila pneumoniae имеет 1112.
M = numel(seqtrachomatis.CDS) N = numel(seqpneumoniae.CDS)
M =
895
N =
1112
Большинство записей CDS содержат трансляцию в аминокислотные последовательности. Первая запись CDS в данных Chlamydia trachomatis - гипотетический белок длиной 591 остаток.
seqtrachomatis.CDS(1)
ans =
struct with fields:
location: 'join(1041920..1042519,1..1176)'
gene: []
product: 'hypothetical protein'
codon_start: '1'
indices: [1041920 1042519 1 1176]
protein_id: 'NP_219502.1'
db_xref: 'GeneID:884145'
note: []
translation: 'MSIRGVGGNGNSRIPSHNGDGSNRRSQNTKGNNKVEDRVCSLYSSRSNENRESPYAVVDVSSMIESTPTSGETTRASRGVFSRFQRGLVRVADKVRRAVQCAWSSVSTRRSSATRAAESGSSSRTARGASSGYREYSPSAARGLRLMFTDFWRTRVLRQTSPMAGVFGNLDVNEARLMAAYTSECADHLEANKLAGPDGVAAAREIAKRWEQRVRDLQDKGAARKLLNDPLGRRTPNYQSKNPGEYTVGNSMFYDGPQVANLQNVDTGFWLDMSNLSDVVLSREIQTGLRARATLEESMPMLENLEERFRRLQETCDAARTEIEESGWTRESASRMEGDEAQGPSRAQQAFQSFVNECNSIEFSFGSFGEHVRVLCARVSRGLAAAGEAIRRCFSCCKGSTHRYAPRDDLSPEGASLAETLARFADDMGIERGADGTYDIPLVDDWRRGVPSIEGEGSDSIYEIMMPIYEVMDMDLETRRSFAVQQGHYQDPRASDYDLPRASDYDLPRSPYPTPPLPPRYQLQNMDVEAGFREAVYASFVAGMYNYVVTQPQERIPNSQQVEGILRDMLTNGSQTFRDLMRRWNREVDRE'
text: [19x58 char]
Четвертая запись CDS относится к гену gatA, который имеет субъединицу продукта глутамил-тРНК амидотрансферазы A. Длина последовательности продукта составляет 491 остаток.
seqtrachomatis.CDS(4)
ans =
struct with fields:
location: '2108..3583'
gene: 'gatA'
product: [2x47 char]
codon_start: '1'
indices: [2108 3583]
protein_id: 'NP_219505.1'
db_xref: 'GeneID:884087'
note: [7x58 char]
translation: 'MYRKSALELRDAVVNRELSVTAITEYFYHRIESHDEQIGAFLSLCKERALLRASRIDDKLAKGDPIGLLAGIPIGVKDNIHITGVKTTCASKMLENFVAPFDSTVVRRIEMEDGILLGKLNMDEFAMGSTTRYSAFHPTNNPWDLERVPGGSSGGSAAAVSARFCPIALGSDTGGSIRQPAAFCGVVGFKPSYGAVSRYGLVAFGSSLDQIGPLTTVVEDVALAMDAFAGRDPKDSTTRDFFKGTFSQALSLEVPKLIGVPRGFLDGLQEDCKENFFEALAVMEREGSRIIDVDLSVLKHAVPVYYIVASAEAATNLARFDGVRYGHRCAQADNMHEMYARSRKEGFGKEVTRRILLGNYVLSAERQNIFYKKGMAVRARLIDAFQAAFERCDVIAMPVCATPAIRDQDVLDPVSLYLQDVYTVAVNLAYLPAISVPSGLSKEGLPLGVQFIGERGSDQQICQVGYSFQEHSQIKQLYPKAVNGLFDGGIE'
text: [26x58 char]
Некоторые из CDS Chlamydophila pneumoniae имеют пустые переводы. Заполните их следующим образом. Сначала найдите все пустые переводы, а затем просмотрите первый пустой перевод.
missingPn = find(cellfun(@isempty,{seqpneumoniae.CDS.translation}));
seqpneumoniae.CDS(missingPn(1))
ans =
struct with fields:
location: 'complement(73364..73477)'
gene: []
product: 'hypothetical protein'
codon_start: '1'
indices: [73477 73364]
protein_id: 'NP_444613.1'
db_xref: 'GeneID:963699'
note: 'hypothetical protein; identified by Glimmer2'
translation: []
text: [10x52 char]
Функция featureparse извлекает элементы, такие как CDS, из структуры последовательности. Затем можно использовать cellfun подавать nt2aa в последовательности с отсутствующими трансляциями.
allCDS = featureparse(seqpneumoniae,'Feature','CDS','Sequence',true); missingSeqs = cellfun(@nt2aa,{allCDS(missingPn).Sequence},'uniform',false); [seqpneumoniae.CDS(missingPn).translation] = deal(missingSeqs{:}); seqpneumoniae.CDS(missingPn(1))
ans =
struct with fields:
location: 'complement(73364..73477)'
gene: []
product: 'hypothetical protein'
codon_start: '1'
indices: [73477 73364]
protein_id: 'NP_444613.1'
db_xref: 'GeneID:963699'
note: 'hypothetical protein; identified by Glimmer2'
translation: 'MLTDQRKHIQMLHKHNSIEIFLSNMVVEVKLFFKTLK*'
text: [10x52 char]
Чтобы сравнить ген gatA в Chlamydia trachomatis со всеми генами CDS в Chlamydophila pneumoniae, поставьте for закольцовывание вокруг nwalign функция. В качестве альтернативы можно использовать локальное выравнивание (swalign).
tic gatAScores = zeros(1,N); for inner = 1:N gatAScores(inner) = nwalign(seqtrachomatis.CDS(4).translation,... seqpneumoniae.CDS(inner).translation); end toc % |tic| and |toc| are used to report how long the calculation takes.
Elapsed time is 1.898762 seconds.
Гистограмма баллов показывает большое количество отрицательных баллов и один очень высокий положительный балл.
hist(gatAScores,100) title(sprintf(['Alignment Scores for Chlamydia trachomatis %s\n',... 'with all CDS in Chlamydophila pneumoniae'],seqtrachomatis.CDS(4).gene))
Как и ожидалось, высокая оценка соответствует гену gatA в Chlamydophila pneumoniae.
[gatABest, gatABestIdx] = max(gatAScores); seqpneumoniae.CDS(gatABestIdx)
ans =
struct with fields:
location: 'complement(838828..840306)'
gene: 'gatA'
product: [2x47 char]
codon_start: '1'
indices: [840306 838828]
protein_id: 'NP_445311.1'
db_xref: 'GeneID:963139'
note: [7x58 char]
translation: 'MYRYSALELAKAVTLGELTATGVTQHFFHRIEEAEGQVGAFISLCKEQALEQAELIDKKRSRGEPLGKLAGVPVGIKDNIHVTGLKTTCASRVLENYQPPFDATVVERIKKEDGIILGKLNMDEFAMGSTTLYSAFHPTHNPWDLSRVPGGSSGGSAAAVSARFCPVALGSDTGGSIRQPAAFCGVVGFKPSYGAVSRYGLVAFASSLDQIGPLANTVEDVALMMDVFSGRDPKDATSREFFRDSFMSKLSTEVPKVIGVPRTFLEGLRDDIRENFFSSLAIFEGEGTHLVDVELDILSHAVSIYYILASAEAATNLARFDGVRYGYRSPQAHTISQLYDLSRGEGFGKEVMRRILLGNYVLSAERQNVYYKKATAVRAKIVKAFRTAFEKCEILAMPVCSSPAFEIGEILDPVTLYLQDIYTVAMNLAYLPAIAVPSGFSKEGLPLGLQIIGQQGQDQQVCQVGYSFQEHAQIKQLFSKRYAKSVVLGGQS'
text: [26x58 char]
Парное выравнивание одного гена Chlamydia trachomatis со всеми генами Chlamydophila pneumoniae занимает чуть менее минуты на машине Intel ® Pentium 4, 2,0 ГГц под управлением Windows ® XP. Для выполнения этого расчета для всех 895 CDS в Chlamydia trachomatis потребуется около 12 часов на одной машине. Раскомментируйте следующий код, если требуется выполнить весь расчет.
scores = zeros(M,N); parfor outer = 1:M theScore = zeros(1,outer); theSeq = seqtrachomatis.CDS(outer).translation; for inner = 1:N theScore(inner) = ... nwalign(theSeq,... seqpneumoniae.CDS(inner).translation); end scores(outer,:) = theScore; end
Обратите внимание на команду parfor используется во внешнем контуре. Если компьютер настроен на выполнение нескольких лабораторий, то внешний цикл будет выполняться параллельно. Для получения полного представления об этой конструкции см. doc parfor.
Распределение баллов для нескольких генов показывает закономерность. CDS (3) Chlamydia trachomatis является геном gatC. Это имеет относительно короткий продукт, аспартил/глутамил-тРНК амидотрансферазную субъединицу С, с только 100 остатками.
gatCScores = zeros(1,N); for inner = 1:N gatCScores(inner) = nwalign(seqtrachomatis.CDS(3).translation,... seqpneumoniae.CDS(inner).translation); end figure hist(gatCScores,100) title(sprintf(['Alignment Scores for Chlamydia trachomatis %s\n',... 'with all CDS in Chlamydophila pneumoniae'],seqtrachomatis.CDS(3).gene)) xlabel('Score');ylabel('Number of genes');
Лучший балл снова соответствует тому же гену в Chlamydophila pneumoniae.
[gatCBest, gatCBestIdx] = max(gatCScores); seqpneumoniae.CDS(gatCBestIdx).product
ans =
2x47 char array
'aspartyl/glutamyl-tRNA amidotransferase subunit'
'C '
CDS (339) Chlamydia trachomatis является геном uvrA. Он имеет очень длинный продукт, субъединицу ABC A, длиной 1786.
uvrAScores = zeros(1,N); for inner = 1:N uvrAScores(inner) = nwalign(seqtrachomatis.CDS(339).translation,... seqpneumoniae.CDS(inner).translation); end figure hist(uvrAScores,100) title(sprintf(['Alignment Scores for Chlamydia trachomatis %s\n',... 'with all CDS in Chlamydophila pneumoniae'],seqtrachomatis.CDS(339).gene)) xlabel('Score');ylabel('Number of genes'); [uvrABest, uvrABestIdx] = max(uvrAScores); seqpneumoniae.CDS(uvrABestIdx)
ans =
struct with fields:
location: '716887..722367'
gene: []
product: 'excinuclease ABC subunit A'
codon_start: '1'
indices: [716887 722367]
protein_id: 'NP_445220.1'
db_xref: 'GeneID:963214'
note: [6x58 char]
translation: 'MKSLPVYVSGIKVRNLKNVSIHFNSEEIVLLTGVSGSGKSSIAFDTLYAAGRKRYISTLPTFFATTITTLPNPKVEEIHGLSPTIAIKQNHFSHYSHATVGSTTELFSHLALLFTLEGQARDPKTKEVLDLYSKEKVLSTIMELSEGVQISILAPLLRKDIAAIHEYAQQGFTKVRCNGTIHPIYSFLTSGIPEDCSVDIVIDTLIKSENNIARLKVSLFTALEFGEGHCSVLSDEELMTFSTKQQIDDVTYTPLTQQLFSPHALESRCSLCQGSGIFISIDNPLLIDENLSIKENCCSFAGNCSSYLYHTIYQALADALNFNLETPWKDLSPEIQNIFLRGKNNLVLPVRLFDQTLGKKNLTYKVWRGVLNDIGDKVRYTTKPSRYLSKGMSAHSCSLCKGTGLGDYASVATWEGKTFTEFQQMSLNNWHVFFSKVKSPSLSIQEILQGLKQRLSFLIDLGLGYLTPNRALATLSGGEQERTAIAKHLGGELFGITYILDEPSIGLHPQDTEKLIGVIKKLRDQGNTVILVEHEERMISLADRIIDIGPGAGIFGGEVLFNGKPEDFLMNSSSLTAKYLRQELTIPIPESREAPTSWLLLTEATIHNLKNLSIRLPLARLIGVTGVSGSGKSSLINNTLVPAIESFLKQENPKNLHFEWGCIGRLIHITRDLPGRSQRSIPLTYIKAFDDIRELFASQPRSLRQGLTKAHFSFNQPQGACIQCQGLGTMTISDDDTPIPCSECQGKRYHSEVLEILYEGKNIADILDMTAYEAEKFFISHPKIHEKIHALCSLRLDYLPLGRPLSTLSGGEIQRLKLAHELLFASPKQTLYVLDEPTTGLHTHDIQALIEVLLSLTYLGHTVLVIEHNMHVVKVCDYVLELGPEGGDLGGYLLASCTPKDLIQLNTPTAKALAPYIEGSLDIPVVKSEPPSSPKSCDILIKDAYQNNLKHIDLALPRNSLIAIAGPGASGKHSLVFDILYASGNIAYAELFPPYIRQGLLKETPLPSVGEVKGLSPVISVRKCSSSNRSYHTIASALGLSNGLEKLFAILGEPFSPLTEEKLSKTTPQTIIDSLLKSYKDDYVTITSPIPLGSDLEIFLQEKQKEGFIKLYSEGNLYDLDERLPLNLIEPAIVIQHTKVSPKNSSSLLSAISVAFSLSSEIWIYISQKKQRKLSYSLGWKDKKGRLYPEITHQLLSSDHPEGRCLTCGGRGEILKISLEEHKEKIAHYTPLEFFSLFFPKSYMKPVQKLLKDENASQPLKLLTTKEFLNFCRGSSEFPGMNALLMEQLDTESDSPLIKPLLALTSCPACKGSGLNDYANYVRINNTSLLDIYQEDATFLESFLNTIGTDDTRSIIQDLMNRLTFISKVGLSYITLGQRQDTLSDGENYRLHLAKKISINLTNIVYLFEEPLSGLHPQDLPTIVQLLKELVANNNTVIATDRSCSLIPHADHAIFLGPGSGPQGGFLMDSDTEVCPSVDLHANVPQTEVCPKAPLSISKANHTRGSDRTLKVNLSIHHIQNLKVSAPLHALVAIGGVSGSGKTSLLLEGFKKQAELLIAKGTTTFSDLVVIDSHPIASSQRSDISTYFDIAPSLRAFYASLTQAKALNISSTMFSTNTKQGQCSDCQGLGYQWIDRAFYALEKRPCPTCSGFRIQPLAQEVLYEGKHFGELLHTPIETVALRFPFIKKIQKPLKALLDIGLGYLPIGQKLSSLSVSEKTALKTAYFLYQTPETPTLFLIDELFSSLDPIKKQHLPEKLRSLINSGHSVIYIDHDVKLLKSADYLIEIGPGSGKQGGKLLFSGSPKDIYASKDSLLKKYICNEELDS'
text: [46x58 char]
На распределение баллов влияет длина последовательностей, причем очень длинные последовательности потенциально имеют гораздо более высокие или более низкие оценки, чем более короткие последовательности. Нормализовать для этого можно несколькими способами. Одним из способов является деление на длину последовательностей.
lnormgatABest = gatABest./length(seqtrachomatis.CDS(4).product) lnormgatCBest = gatCBest./length(seqtrachomatis.CDS(3).product) lnormuvrABest = uvrABest./length(seqtrachomatis.CDS(339).product)
lnormgatABest =
16.8794
lnormgatCBest =
2.2695
lnormuvrABest =
78.9615
Альтернативным методом нормализации является использование показателя самосогласования, то есть балла от выравнивания последовательности с самим собой.
gatASelf = nwalign(seqtrachomatis.CDS(4).translation,... seqtrachomatis.CDS(4).translation); gatCSelf = nwalign(seqtrachomatis.CDS(3).translation,... seqtrachomatis.CDS(3).translation); uvrASelf = nwalign(seqtrachomatis.CDS(339).translation,... seqtrachomatis.CDS(339).translation); normgatABest = gatABest./gatASelf normgatCBest = gatCBest./gatCSelf normuvrABest = uvrABest./uvrASelf
normgatABest =
0.7380
normgatCBest =
0.5212
normuvrABest =
0.5253
Расчет выравнивания «все против всех» не только занимает много времени, но и создает большую матрицу баллов. Если вы ищете похожие гены для разных видов, то интересные оценки являются положительными оценками, которые указывают на хорошее выравнивание. Однако большинство из этих баллов являются отрицательными, и фактические значения не особенно полезны для этого типа исследований. Разреженные матрицы позволяют хранить интересные значения более эффективным способом.
Разреженная матрица, spScores, в MAT-файле chlamydia.mat содержит положительные значения от всех по отношению ко всем вычислениям парной центровки, нормализованным по баллу самонастройки.
load('chlamydia.mat','spScores')
С помощью матрицы баллов можно посмотреть на распределение баллов генов Chlamydophila pneumoniae, выровненных с Chlamydia trachomatis, и на обратное это, гены Chlamydia trachomatis, выровненные с генами Chlamydophila pneumoniae
figure subplot(2,1,1) hist(max(spScores),100) title('Highest Alignment Scores for Chlamydophila pneumoniae Genes') xlabel('Score');ylabel('Number of genes'); subplot(2,1,2) hist(max(spScores,[],2),100) title('Highest Alignment Scores for Chlamydia trachomatis Genes') xlabel('Score');ylabel('Number of genes');
Помните, что существует 1112 CDS в Chlamydophila pneumoniae и только 895 в Chlamydia trachomatis. Большое количество нулевых баллов в верхней гистограмме указывает на то, что многие из дополнительных CDS у Chlamydophila pneumoniae не имеют хороших совпадений у Chlamydia trachomatis.
Другим способом визуализации данных является просмотр позиций точек в матрице оценок, которые являются положительными. Разреженная функция spy является простым способом быстрого просмотра dotplots матриц. Это показывает какую-то интересную структуру на позициях высоких голевых матчей.
figure
spy(spScores > 0)
title(sprintf('Dot Plot of High-Scoring Alignments.\nNormalized Threshold = 0'))
Поднимите порог чуть выше, чтобы увидеть на сюжете четкие диагональные линии.
spy(spScores >.1)
title(sprintf('Dot Plot of High-Scoring Alignments.\nNormalized Threshold = 0.1'))
Помните, что это циркулярные геномы, и кажется, что отправные точки в GenBank произвольны. Переставьте матрицу оценок так, чтобы лучшее совпадение первого CDS у Chlamydophila pneumoniae было в первом ряду, чтобы увидеть четкий диагональный график. Это показывает синтению между двумя организмами.
[bestScore bestMatch] = max(spScores(:,1));
spy(spScores([bestMatch:end 1:bestMatch-1],:)>.1);
title('Synteny Plot of Chlamydophila pneumoniae and Chlamydia trachomatis')
Гены в различных геномах, которые связаны друг с другом, считаются гомологичными. Сходство может быть по видообразованию (ортологичные гены) или по репликации (паралогозные гены). Наличие матрицы оценки позволяет искать оба типа отношений.
Наиболее очевидным способом найти ортологи является поиск наивысшей оценки пар для каждого гена. Если оценка значительна, то эти лучшие обратные пары с большой вероятностью будут ортологичными.
[bestScores, bestIndices] = max(spScores);
Переменная bestIndices содержит индекс лучших обратных пар для генов в Chlamydophila pneumoniae. Сортируйте лучшие баллы и создайте таблицу, чтобы сравнить описание лучших обратных пар и обнаружить очень высокое сходство между наивысшими баллами лучших обратных пар.
[orderedScores, permScores] = sort(full(bestScores),'descend'); matches = [num2cell(orderedScores)',num2cell(bestIndices(permScores))',... num2cell((permScores))',... {seqtrachomatis.CDS(bestIndices(permScores)).product;... seqpneumoniae.CDS((permScores)).product; }']; for count = 1:7 fprintf(['Score %f\nChlamydia trachomatis Gene : %s\n',... 'Chlamydophila pneumoniae Gene : %s\n\n'],... matches{count,1}, matches{count,4}, matches{count,5}) end
Score 0.982993 Chlamydia trachomatis Gene : 50S ribosomal protein L36 Chlamydophila pneumoniae Gene : 50S ribosomal protein L36 Score 0.981818 Chlamydia trachomatis Gene : 30S ribosomal protein S15 Chlamydophila pneumoniae Gene : 30S ribosomal protein S15 Score 0.975422 Chlamydia trachomatis Gene : integration host factor alpha-subunit Chlamydophila pneumoniae Gene : integration host factor beta-subunit Score 0.971647 Chlamydia trachomatis Gene : 50S ribosomal protein L16 Chlamydophila pneumoniae Gene : 50S ribosomal protein L16 Score 0.970105 Chlamydia trachomatis Gene : 30S ribosomal protein S10 Chlamydophila pneumoniae Gene : 30S ribosomal protein S10 Score 0.969554 Chlamydia trachomatis Gene : rod shape-determining protein MreB Chlamydophila pneumoniae Gene : rod shape-determining protein MreB Score 0.953654 Chlamydia trachomatis Gene : hypothetical protein Chlamydophila pneumoniae Gene : hypothetical protein
Редактор переменных можно использовать для просмотра данных в формате электронной таблицы.
open('matches')
Сравните описания, чтобы увидеть, что большинство лучших взаимных пар имеют идентичные описания.
exactMatches = strcmpi(matches(:,4),matches(:,5)); sum(exactMatches)
ans = 808
Возможно, более интересны лучшие обратные пары, где описания не идентичны. Некоторые из них являются просто различиями в том, как описывается один и тот же ген, но другие показывают совершенно разные описания.
mismatches = matches(~exactMatches,:); for count = 1:7 fprintf(['Score %f\nChlamydia trachomatis Gene : %s\n',... 'Chlamydophila pneumoniae Gene : %s\n\n'],... mismatches{count,1}, mismatches{count,4}, mismatches{count,5}) end
Score 0.975422 Chlamydia trachomatis Gene : integration host factor alpha-subunit Chlamydophila pneumoniae Gene : integration host factor beta-subunit Score 0.929565 Chlamydia trachomatis Gene : low calcium response D Chlamydophila pneumoniae Gene : type III secretion inner membrane protein SctV Score 0.905000 Chlamydia trachomatis Gene : NrdR family transcriptional regulator Chlamydophila pneumoniae Gene : transcriptional regulator NrdR Score 0.903226 Chlamydia trachomatis Gene : Yop proteins translocation protein S Chlamydophila pneumoniae Gene : type III secretion inner membrane protein SctS Score 0.896212 Chlamydia trachomatis Gene : ATP-dependent protease ATP-binding subunit ClpX Chlamydophila pneumoniae Gene : ATP-dependent protease ATP-binding protein ClpX Score 0.890705 Chlamydia trachomatis Gene : ribonuclease E Chlamydophila pneumoniae Gene : ribonuclease G Score 0.884234 Chlamydia trachomatis Gene : ClpC protease ATPase Chlamydophila pneumoniae Gene : ATP-dependent Clp protease ATP-binding protein
Просмотр данных для несовпадений.
open('mismatches')
Как только у вас есть матрица оценки, это открывает много возможностей для дальнейшего расследования. Например, вы можете искать CDS, где есть несколько высоких баллов, обратных CDS. Дополнительные идеи см. у Кристианини и Хана [1].
[1] Криштианини, Н. и Хан, М.В., «Введение в вычислительную геномику: подход к тематическим исследованиям», Cambridge University Press, 2007.