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

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

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

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

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