exponenta event banner

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

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

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

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

Консенсусная последовательность генома митохондрий человека имеет номер присоединения 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 также полезно для идентификации генов и для оценки температуры отжига последовательности ДНК.

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

Теперь вы будете использовать некоторые функции статистики последовательностей в 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. GeneticCode опции для seqshoworfs функция позволяет снова взглянуть на ORF, но на этот раз с митохондриальным генетическим кодом позвоночных.

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

Обратите внимание, что теперь на третьей рамке считывания есть два гораздо больших ORF: один начинается с позиции 4470, а другой начинается с 5904. Они соответствуют генам ND2 (субъединица 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 функция извлечения, когда это возможно, последовательностей ДНК, соответствующих каждому признаку. 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), который позволяет легко создавать графики различных свойств, таких как гидрофобность, белковой последовательности. Выберите меню «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, чем в других. Вы можете проверить, в какие аминокислоты транслируются эти кодоны, используя 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. и Drouin, J., «Различный генетический код в митохондриях человека», Nature, 282 (5735): 189-94, 1979.