Вычисление и визуализация статистики последовательности
В этом примере показано, как использовать основные методы манипуляции с последовательностями и вычисляет некоторую полезную статистику последовательностей. Это также иллюстрирует, как искать кодирующие области (такие как белки) и преследовать последующий анализ из них.
Митохондриальный геном человека
В этом примере вы исследуете последовательность ДНК митохондрий человека. Митохондрии - это структуры, называемые органеллы, которые обнаруживаются в цитоплазме камеры в сотнях-тысячах на каждую камеру. Митохондрии, как правило, являются основным центром производства энергии в эукариотах, они помогают деградировать жиры и сахара.
Консенсусная последовательность генома митохондрий человека имеет число NC_012920 присоединения. Можно getgenbank функция для получения последней аннотированной последовательности из GenBank ® в рабочую область MATLAB ®.
mitochondria_gbk = getgenbank('NC_012920');
Для вашего удобства ранее загруженная последовательность включается в MAT-файл. Обратите внимание, что данные в общедоступных хранилищах часто кураторятся и обновляются; поэтому результаты этого примера могут несколько отличаться при использовании современных наборов данных.
load mitochondria
Скопируйте только последовательность ДНК в новую переменную mitochondria. Вы можете получить доступ к частям последовательности ДНК с помощью регулярных команд индексации MATLAB.
mitochondria = mitochondria_gbk.Sequence; mitochondria_length = length(mitochondria) first_300_bases = seqdisp(mitochondria(1:300))
mitochondria_length =
16569
first_300_bases =
5×70 char array
' 1 GATCACAGGT CTATCACCCT ATTAACCACT CACGGGAGCT CTCCATGCAT TTGGTATTTT'
' 61 CGTCTGGGGG GTATGCACGC GATAGCATTG CGAGACGCTG GAGCCGGAGC ACCCTATGTC'
'121 GCAGTATCTG TCTTTGATTC CTGCCTCATC CTATTATTTA TCGCACCTAC GTTCAATATT'
'181 ACAGGCGAAC ATACTTACTA AAGTGTGTTA ATTAATTAAT GCTTGTAGGA CATAATAATA'
'241 ACAATTGAAT GTCTGCACAG CCACTTTCCA CACAGACATC ATAACAAAAA ATTTCCACCA'
Можно посмотреть на состав нуклеотидов с ntdensity функция.
figure ntdensity(mitochondria)
Это показывает, что геном митохондрий богат A-T. Содержание GC иногда используется для классификации организмов по систематике, оно может варьироваться между различными видами от ~ 30% до ~ 70%. Измерение содержимого GC также полезно для идентификации генов и для оценки температуры отжига последовательности ДНК.
Вычисление статистики последовательности
Теперь вы будете использовать некоторые функции статистики последовательностей в Bioinformatics Toolbox™, чтобы посмотреть на различные свойства митохондриального генома человека. Количество основ всей последовательности можно подсчитать с помощью basecount функция.
bases = basecount(mitochondria)
bases =
struct with fields:
A: 5124
C: 5181
G: 2169
T: 4094
Они находятся на цепи 5 '-3'. Вы можете посмотреть на обратный случай дополнения, используя seqrcomplement функция.
compBases = basecount(seqrcomplement(mitochondria))
compBases =
struct with fields:
A: 4094
C: 2169
G: 5181
T: 5124
Как ожидалось, отсчеты основы на обратной цепи дополнения комплементарны отсчетам на цепи 5 '-3'.
Можно использовать опцию графика, чтобы basecount отображение круговой диаграммы распределения основ.
figure basecount(mitochondria,'chart','pie'); title('Distribution of Nucleotide Bases for Human Mitochondrial Genome');
Теперь посмотрите на димеры в последовательности и отобразите информацию на столбчатой диаграмме с помощью dimercount.
figure dimers = dimercount(mitochondria,'chart','bar') title('Mitochondrial Genome Dimer Histogram');
dimers =
struct with fields:
AA: 1604
AC: 1495
AG: 795
AT: 1230
CA: 1534
CC: 1771
CG: 435
CT: 1440
GA: 613
GC: 711
GG: 425
GT: 419
TA: 1373
TC: 1204
TG: 513
TT: 1004
Исследование открытых систем координат чтения (ORF)
В нуклеотидной последовательности очевидная вещь, которую нужно искать, есть ли какие-либо открытые системы координат считывания. ORF представляет собой любую последовательность ДНК или РНК, которая потенциально может быть переведена в белок. Функция seqshoworfs может использоваться для визуализации ORF в последовательности.
Примечание: В руководстве HTML показана только первая страница выхода, однако при запуске примера вы сможете просмотреть полный геном митохондрий с помощью полосы прокрутки на рисунке.
seqshoworfs(mitochondria);
Если сравнить этот выход с генами, показанными на странице NCBI, то, по-видимому, немного меньше ORF, и, следовательно, меньше генов, чем ожидалось.
Митохондрии позвоночных не используют Стандартный генетический код, поэтому некоторые кодоны имеют разное значение в митохондриальных геномах. Для получения дополнительной информации об использовании различных генетических кодов в MATLAB смотрите справку по функции geneticcode. The GeneticCode опция для seqshoworfs функция позволяет вам снова посмотреть на ORF, но на этот раз с митохондриальным генетическим кодом позвоночных.
В последовательности митохондриальной ДНК человека некоторые гены также начинаются альтернативными стартовыми кодонами [1]. Используйте AlternativeStartCodons опция для seqshoworfs функция для поиска также этих ORF.
Заметьте, что теперь на третьей системе координат считывания есть два намного больших ORF: Один, начиная с позиции 4470, и другой, начиная с 5904. Они соответствуют генам ND2 (NADH-дегидрогеназа подмодуля 2) и COX1 (цитохром с оксидазой подмодуля I).
orfs = seqshoworfs(mitochondria,'GeneticCode','Vertebrate Mitochondrial',... 'AlternativeStartCodons',true)
orfs =
1×3 struct array with fields:
Start
Stop
Проверка аннотированных функций
Можно также посмотреть на все функции, которые были аннотированы к геному митохондрий человека. Исследуйте полный ввод GenBank mitochondria_gbk с featureparse функция. В частности, можно исследовать аннотированные последовательности кодирования (CDS) и сравнить их с ранее найденными ORF. Используйте Sequence опция для featureparse функция для извлечения, когда это возможно, последовательностей ДНК, соответствующих каждому признаку. The featureparse функция будет дополнять части исходной последовательности, когда это необходимо.
features = featureparse(mitochondria_gbk,'Sequence',true) coding_sequences = features.CDS; coding_sequences_id = sprintf('%s ',coding_sequences.gene)
features =
struct with fields:
source: [1×1 struct]
D_loop: [1×1 struct]
gene: [1×37 struct]
tRNA: [1×22 struct]
rRNA: [1×2 struct]
STS: [1×28 struct]
misc_feature: [1×1 struct]
CDS: [1×13 struct]
coding_sequences_id =
'ND1 ND2 COX1 COX2 ATP8 ATP6 COX3 ND3 ND4L ND4 ND5 ND6 CYTB '
ND2CDS = coding_sequences(2) % ND2 is in the 2nd position COX1CDS = coding_sequences(3) % COX1 is in the 3rd position
ND2CDS =
struct with fields:
Location: '4470..5511'
Indices: [4470 5511]
gene: 'ND2'
gene_synonym: 'MTND2'
note: 'TAA stop codon is completed by the addition of 3' A residues to the mRNA'
codon_start: '1'
transl_except: '(pos:5511,aa:TERM)'
transl_table: '2'
product: 'NADH dehydrogenase subunit 2'
protein_id: 'YP_003024027.1'
db_xref: {'GI:251831108' 'GeneID:4536' 'HGNC:7456' 'MIM:516001'}
translation: 'MNPLAQPVIYSTIFAGTLITALSSHWFFTWVGLEMNMLAFIPVLTKKMNPRSTEAAIKYFLTQATASMILLMAILFNNMLSGQWTMTNTTNQYSSLMIMMAMAMKLGMAPFHFWVPEVTQGTPLTSGLLLLTWQKLAPISIMYQISPSLNVSLLLTLSILSIMAGSWGGLNQTQLRKILAYSSITHMGWMMAVLPYNPNMTILNLTIYIILTTTAFLLLNLNSSTTTLLLSRTWNKLTWLTPLIPSTLLSLGGLPPLTGFLPKWAIIEEFTKNNSLIIPTIMATITLLNLYFYLRLIYSTSITLLPMSNNVKMKWQFEHTKPTPFLPTLIALTTLLLPISPFMLMIL'
Sequence: 'attaatcccctggcccaacccgtcatctactctaccatctttgcaggcacactcatcacagcgctaagctcgcactgattttttacctgagtaggcctagaaataaacatgctagcttttattccagttctaaccaaaaaaataaaccctcgttccacagaagctgccatcaagtatttcctcacgcaagcaaccgcatccataatccttctaatagctatcctcttcaacaatatactctccggacaatgaaccataaccaatactaccaatcaatactcatcattaataatcataatagctatagcaataaaactaggaatagccccctttcacttctgagtcccagaggttacccaaggcacccctctgacatccggcctgcttcttctcacatgacaaaaactagcccccatctcaatcatataccaaatctctccctcactaaacgtaagccttctcctcactctctcaatcttatccatcatagcaggcagttgaggtggattaaaccaaacccagctacgcaaaatcttagcatactcctcaattacccacataggatgaataatagcagttctaccgtacaaccctaacataaccattcttaatttaactatttatattatcctaactactaccgcattcctactactcaacttaaactccagcaccacgaccctactactatctcgcacctgaaacaagctaacatgactaacacccttaattccatccaccctcctctccctaggaggcctgcccccgctaaccggctttttgcccaaatgggccattatcgaagaattcacaaaaaacaatagcctcatcatccccaccatcatagccaccatcaccctccttaacctctacttctacctacgcctaatctactccacctcaatcacactactccccatatctaacaacgtaaaaataaaatgacagtttgaacatacaaaacccaccccattcctccccacactcatcgcccttaccacgctactcctacctatctccccttttatactaataatcttat'
COX1CDS =
struct with fields:
Location: '5904..7445'
Indices: [5904 7445]
gene: 'COX1'
gene_synonym: 'COI; MTCO1'
note: 'cytochrome c oxidase I'
codon_start: '1'
transl_except: []
transl_table: '2'
product: 'cytochrome c oxidase subunit I'
protein_id: 'YP_003024028.1'
db_xref: {'GI:251831109' 'GeneID:4512' 'HGNC:7419' 'MIM:516030'}
translation: 'MFADRWLFSTNHKDIGTLYLLFGAWAGVLGTALSLLIRAELGQPGNLLGNDHIYNVIVTAHAFVMIFFMVMPIMIGGFGNWLVPLMIGAPDMAFPRMNNMSFWLLPPSLLLLLASAMVEAGAGTGWTVYPPLAGNYSHPGASVDLTIFSLHLAGVSSILGAINFITTIINMKPPAMTQYQTPLFVWSVLITAVLLLLSLPVLAAGITMLLTDRNLNTTFFDPAGGGDPILYQHLFWFFGHPEVYILILPGFGMISHIVTYYSGKKEPFGYMGMVWAMMSIGFLGFIVWAHHMFTVGMDVDTRAYFTSATMIIAIPTGVKVFSWLATLHGSNMKWSAAVLWALGFIFLFTVGGLTGIVLANSSLDIVLHDTYYVVAHFHYVLSMGAVFAIMGGFIHWFPLFSGYTLDQTYAKIHFTIMFIGVNLTFFPQHFLGLSGMPRRYSDYPDAYTTWNILSSVGSFISLTAVMLMIFMIWEAFASKRKVLMVEEPSMNLEWLYGCPPPYHTFEEPVYMKS'
Sequence: 'atgttcgccgaccgttgactattctctacaaaccacaaagacattggaacactatacctattattcggcgcatgagctggagtcctaggcacagctctaagcctccttattcgagccgagctgggccagccaggcaaccttctaggtaacgaccacatctacaacgttatcgtcacagcccatgcatttgtaataatcttcttcatagtaatacccatcataatcggaggctttggcaactgactagttcccctaataatcggtgcccccgatatggcgtttccccgcataaacaacataagcttctgactcttacctccctctctcctactcctgctcgcatctgctatagtggaggccggagcaggaacaggttgaacagtctaccctcccttagcagggaactactcccaccctggagcctccgtagacctaaccatcttctccttacacctagcaggtgtctcctctatcttaggggccatcaatttcatcacaacaattatcaatataaaaccccctgccataacccaataccaaacgcccctcttcgtctgatccgtcctaatcacagcagtcctacttctcctatctctcccagtcctagctgctggcatcactatactactaacagaccgcaacctcaacaccaccttcttcgaccccgccggaggaggagaccccattctataccaacacctattctgatttttcggtcaccctgaagtttatattcttatcctaccaggcttcggaataatctcccatattgtaacttactactccggaaaaaaagaaccatttggatacataggtatggtctgagctatgatatcaattggcttcctagggtttatcgtgtgagcacaccatatatttacagtaggaatagacgtagacacacgagcatatttcacctccgctaccataatcatcgctatccccaccggcgtcaaagtatttagctgactcgccacactccacggaagcaatatgaaatgatctgctgcagtgctctgagccctaggattcatctttcttttcaccgtaggtggcctgactggcattgtattagcaaactcatcactagacatcgtactacacgacacgtactacgttgtagcccacttccactatgtcctatcaataggagctgtatttgccatcataggaggcttcattcactgatttcccctattctcaggctacaccctagaccaaacctacgccaaaatccatttcactatcatattcatcggcgtaaatctaactttcttcccacaacactttctcggcctatccggaatgccccgacgttactcggactaccccgatgcatacaccacatgaaacatcctatcatctgtaggctcattcatttctctaacagcagtaatattaataattttcatgatttgagaagccttcgcttcgaagcgaaaagtcctaatagtagaagaaccctccataaacctggagtgactatatggatgccccccaccctaccacacattcgaagaacccgtatacataaaatctaga'
Создайте карту с указанием всех функций, обнаруженных в этой записи GenBank, используя featureview функция.
[h,l] = featureview(mitochondria_gbk,{'CDS','tRNA','rRNA','D_loop'},...
[2 1 2 2 2],'Fontsize',9);
legend(h,l,'interpreter','none');
title('Homo sapiens mitochondrion, complete genome')
Извлечение и анализ ND2 и COX1 белков
Можно транслировать последовательности ДНК, которые кодируют белки ND2 и COX1 при помощи nt2aa функция. Снова, GeneticCode опция должна использоваться для определения генетического кода митохондрий позвоночных.
ND2 = nt2aa(ND2CDS,'GeneticCode','Vertebrate Mitochondrial'); disp(seqdisp(ND2))
1 MNPLAQPVIY STIFAGTLIT ALSSHWFFTW VGLEMNMLAF IPVLTKKMNP RSTEAAIKYF 61 LTQATASMIL LMAILFNNML SGQWTMTNTT NQYSSLMIMM AMAMKLGMAP FHFWVPEVTQ 121 GTPLTSGLLL LTWQKLAPIS IMYQISPSLN VSLLLTLSIL SIMAGSWGGL NQTQLRKILA 181 YSSITHMGWM MAVLPYNPNM TILNLTIYII LTTTAFLLLN LNSSTTTLLL SRTWNKLTWL 241 TPLIPSTLLS LGGLPPLTGF LPKWAIIEEF TKNNSLIIPT IMATITLLNL YFYLRLIYST 301 SITLLPMSNN VKMKWQFEHT KPTPFLPTLI ALTTLLLPIS PFMLMIL
COX1 = nt2aa(COX1CDS,'GeneticCode','Vertebrate Mitochondrial'); disp(seqdisp(COX1))
1 MFADRWLFST NHKDIGTLYL LFGAWAGVLG TALSLLIRAE LGQPGNLLGN DHIYNVIVTA 61 HAFVMIFFMV MPIMIGGFGN WLVPLMIGAP DMAFPRMNNM SFWLLPPSLL LLLASAMVEA 121 GAGTGWTVYP PLAGNYSHPG ASVDLTIFSL HLAGVSSILG AINFITTIIN MKPPAMTQYQ 181 TPLFVWSVLI TAVLLLLSLP VLAAGITMLL TDRNLNTTFF DPAGGGDPIL YQHLFWFFGH 241 PEVYILILPG FGMISHIVTY YSGKKEPFGY MGMVWAMMSI GFLGFIVWAH HMFTVGMDVD 301 TRAYFTSATM IIAIPTGVKV FSWLATLHGS NMKWSAAVLW ALGFIFLFTV GGLTGIVLAN 361 SSLDIVLHDT YYVVAHFHYV LSMGAVFAIM GGFIHWFPLF SGYTLDQTYA KIHFTIMFIG 421 VNLTFFPQHF LGLSGMPRRY SDYPDAYTTW NILSSVGSFI SLTAVMLMIF MIWEAFASKR 481 KVLMVEEPSM NLEWLYGCPP PYHTFEEPVY MKS*
Вы можете получить более полное представление об аминокислоте, содержимом с aacount.
figure subplot(2,1,1) ND2aaCount = aacount(ND2,'chart','bar'); title('Histogram of Amino Acid Count for the ND2 Protein'); subplot(2,1,2) COX1aaCount = aacount(COX1,'chart','bar'); title('Histogram of Amino Acid Count for the COX1 Protein');
Заметьте высокое содержимое лейцина, треонина и изолейцина, а также отсутствие цистеина или аспарагиновой кислоты.
Можно использовать atomiccomp и molweight функции для вычисления большего количества свойств белка ND2.
ND2AtomicComp = atomiccomp(ND2) ND2MolWeight = molweight(ND2)
ND2AtomicComp =
struct with fields:
C: 1818
H: 2882
N: 420
O: 471
S: 25
ND2MolWeight =
3.8960e+04
Для дальнейшего исследования свойств ND2 белка используйте proteinplot. Это графический пользовательский интерфейс (GUI), который позволяет легко создавать графики различных свойств, таких как гидрофобность, белковой последовательности. Щелкните в меню «Редактировать», чтобы создать новые свойства, изменить существующие значения свойств или, чтобы настроить параметры сглаживания. Нажмите на меню «Справка» в графическом интерфейсе пользователя для получения дополнительной информации о том, как использовать инструмент.
proteinplot(ND2)
Можно также программно создать графики различных свойств последовательности с помощью proteinpropplot.
figure proteinpropplot(ND2,'PropertyTitle','Parallel beta strand')
Вычисление частоты кодона с использованием всех генов в геноме митохондрий человека
The codoncount функция рассчитывает количество вхождений каждого кодона в последовательности и отображает форматированную таблицу результата.
codoncount(ND2CDS)
AAA - 10 AAC - 14 AAG - 2 AAT - 6 ACA - 11 ACC - 24 ACG - 3 ACT - 5 AGA - 0 AGC - 4 AGG - 0 AGT - 1 ATA - 23 ATC - 24 ATG - 1 ATT - 8 CAA - 8 CAC - 3 CAG - 2 CAT - 1 CCA - 4 CCC - 12 CCG - 2 CCT - 5 CGA - 0 CGC - 3 CGG - 0 CGT - 1 CTA - 26 CTC - 18 CTG - 4 CTT - 7 GAA - 5 GAC - 0 GAG - 1 GAT - 0 GCA - 8 GCC - 7 GCG - 1 GCT - 4 GGA - 5 GGC - 7 GGG - 0 GGT - 1 GTA - 3 GTC - 2 GTG - 0 GTT - 3 TAA - 0 TAC - 8 TAG - 0 TAT - 2 TCA - 7 TCC - 11 TCG - 1 TCT - 4 TGA - 10 TGC - 0 TGG - 1 TGT - 0 TTA - 8 TTC - 7 TTG - 1 TTT - 8
Заметьте, что в гене ND2 больше кодонов CTA, ATC и ACC, чем других. Можно проверить, в какие аминокислоты переводятся эти кодоны, используя nt2aa и aminolookup функций.
CTA_aa = aminolookup('code',nt2aa('CTA')) ATC_aa = aminolookup('code',nt2aa('ATC')) ACC_aa = aminolookup('code',nt2aa('ACC'))
CTA_aa =
'Leu Leucine
'
ATC_aa =
'Ile Isoleucine
'
ACC_aa =
'Thr Threonine
'
Чтобы вычислить частоту кодона для всех генов, вы можете объединить их в одну последовательность перед использованием функции codoncount. Вы должны убедиться, что кодоны полны (три нуклеотида каждый), поэтому считанная система координат последовательности не теряется при конкатенации.
numCDS = numel(coding_sequences); CDS = cell(numCDS,1); for i = 1:numCDS seq = coding_sequences(i).Sequence; CDS{i} = seq(1:3*floor(length(seq)/3)); end allCDS = [CDS{:}]; codoncount(allCDS)
AAA - 85 AAC - 132 AAG - 10 AAT - 32 ACA - 134 ACC - 155 ACG - 10 ACT - 52 AGA - 1 AGC - 39 AGG - 1 AGT - 14 ATA - 167 ATC - 196 ATG - 40 ATT - 124 CAA - 82 CAC - 79 CAG - 8 CAT - 18 CCA - 52 CCC - 119 CCG - 7 CCT - 41 CGA - 28 CGC - 26 CGG - 2 CGT - 7 CTA - 276 CTC - 167 CTG - 45 CTT - 65 GAA - 64 GAC - 51 GAG - 24 GAT - 15 GCA - 80 GCC - 124 GCG - 8 GCT - 43 GGA - 67 GGC - 87 GGG - 34 GGT - 24 GTA - 70 GTC - 48 GTG - 18 GTT - 31 TAA - 3 TAC - 89 TAG - 2 TAT - 46 TCA - 83 TCC - 99 TCG - 7 TCT - 32 TGA - 93 TGC - 17 TGG - 11 TGT - 5 TTA - 73 TTC - 139 TTG - 18 TTT - 77
Используйте figure опция для codoncount функция для отображения тепловой карты с частотой кодона. Используйте geneticcode опция наложения сетки на рисунок, которая группирует синонимические кодоны согласно митохондриальному генетическому коду позвоночных. Наблюдайте конкретное смещение лейцина (кодоны 'CTN').
figure count = codoncount(allCDS,'figure',true,'geneticcode','Vertebrate Mitochondrial'); title('Human Mitochondrial Genome Codon Frequency')
close all
Ссылки
[1] Barrell, B.G., Bankier, A.T. and Drouin, J., «A different genetic кода in human mitochondria», Nature, 282 (5735): 189-94, 1979.