Вычисление и визуализация статистики последовательности

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

Человеческий митохондриальный геном

В этом примере вы исследуете последовательность DNA человеческих митохондрий. Митохондрии являются структурами, названными органоидами, которые найдены в цитоплазме ячейки в сотнях к тысячам для каждой ячейки. Митохондрии обычно являются крупнейшим центром выработки энергии у эукариотов, они помогают ухудшить жиры и сахар.

Последовательность согласия человеческого генома митохондрий имеет инвентарный номер NC_012920. Вы можете getgenbank функционировать, чтобы получить последнюю аннотируемую последовательность от GenBank® в рабочую область MATLAB®.

mitochondria_gbk = getgenbank('NC_012920');

Для вашего удобства ранее загруженная последовательность включена в MAT-файл. Обратите внимание на то, что данные в общедоступных репозиториях часто курируются и обновляются; поэтому результаты этого примера могут немного отличаться, когда вы используете актуальные наборы данных.

load mitochondria

Скопируйте только последовательность DNA в новую переменную mitochondria. Можно получить доступ к частям последовательности DNA при помощи обычного MATLAB, индексирующего команды.

mitochondria = mitochondria_gbk.Sequence;
mitochondria_length = length(mitochondria)
first_300_bases = seqdisp(mitochondria(1:300))
mitochondria_length =

       16569


first_300_bases =

  5x70 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 также полезно для идентификации генов и для оценки температуры отжига последовательности DNA.

Вычисление статистики последовательности

Теперь, вы будете использовать некоторые функции статистики последовательности в 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

Исследование открытых рамок считывания (ORFs)

В последовательности нуклеотида очевидная вещь искать состоит в том, если существуют какие-либо открытые рамки считывания. ORF является любой последовательностью DNA или RNA, который может быть потенциально переведен в белок. Функциональный seqshoworfs может использоваться, чтобы визуализировать ORFs в последовательности.

Примечание: В руководстве по HTML только первую страницу вывода показывают, однако при выполнении примера, вы сможете осмотреть полный митохондриальный геном с помощью полосы прокрутки на фигуре.

seqshoworfs(mitochondria);

Если вы выдерживаете сравнение этот вывод к генам, показанным на странице NCBI, кажется, существует немного меньше ORFs, и следовательно меньше генов, чем ожидалось.

Позвоночные митохондрии не используют Стандартный генетический код, таким образом, некоторые кодоны имеют различное значение в митохондриальных геномах. Для получения дополнительной информации об использовании различных генетических кодов в MATLAB смотрите справку для функционального geneticcode. Опция GeneticCode к функции seqshoworfs позволяет вам смотреть на ORFs снова, но на этот раз с позвоночным митохондриальным генетическим кодом.

В человеческой митохондриальной последовательности DNA некоторые гены также начаты альтернативными кодонами запуска [1]. Используйте опцию AlternativeStartCodons для функции seqshoworfs, чтобы искать также эти ORFs.

Заметьте, что существует теперь два намного больших ORFs на третьей рамке считывания: Один запуск в положении 4470 и другой запуск в 5 904. Они соответствуют ND2 (подблок дегидрогеназы NADH 2) и COX1 (цитохром c подблок оксидазы I) гены.

orfs = seqshoworfs(mitochondria,'GeneticCode','Vertebrate Mitochondrial',...
        'AlternativeStartCodons',true)
orfs = 

  1x3 struct array with fields:

    Start
    Stop

Осмотр аннотируемых функций

Можно также посмотреть на все функции, которые были аннотированы к человеческому митохондриальному геному. Исследуйте полную запись GenBank mitochondria_gbk с функцией featureparse. Особенно, можно исследовать аннотируемые последовательности кодирования (CDS) и сравнить их с ORFs, ранее найденным. Используйте опцию Sequence для функции featureparse, чтобы извлечь, если это возможно, последовательности DNA, соответствующие к каждой функции. Функция featureparse дополнит части исходной последовательности в надлежащих случаях.

features = featureparse(mitochondria_gbk,'Sequence',true)
coding_sequences = features.CDS;
coding_sequences_id = sprintf('%s ',coding_sequences.gene)
features = 

  struct with fields:

          source: [1x1 struct]
          D_loop: [1x1 struct]
            gene: [1x37 struct]
            tRNA: [1x22 struct]
            rRNA: [1x2 struct]
             STS: [1x28 struct]
    misc_feature: [1x1 struct]
             CDS: [1x13 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

Можно перевести последовательности DNA, что код для 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), который позволяет вам легко создавать графики различных свойств, такие как гидрофобность, последовательности белка. Нажмите на меню "Edit", чтобы создать новые свойства, изменить существующие значения свойств, или, настроить параметры сглаживания. Нажмите на меню "Help" в графический интерфейсе пользователя для получения дополнительной информации о том, как использовать инструмент.

proteinplot(ND2)

Можно также программно создать графики различных свойств последовательности с помощью proteinpropplot.

figure
proteinpropplot(ND2,'PropertyTitle','Parallel beta strand')

Вычисление Частоты Кодона с помощью всех Генов в Человеческом Митохондриальном Геноме

Функция 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, чем другие. Можно проверять, какие аминокислоты эти кодоны переводятся в использование функции aminolookup и nt2aa.

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')

Ссылки

[1] Баррель, B.G., Bankier, A.T. и Drouin, J., "Различный генетический код в человеческих митохондриях", Природа, 282 (5735):189-94, 1979.