Exponenta Banner

Оценка значимости выравнивания

Открыть скрипт

Этот пример показывает метод, который может использоваться, чтобы исследовать значимость выравниваний последовательности. Количество тождеств или срабатываний в выравнивании не является четким показателем значительного выравнивания. Сочетание последовательности из трассы будет иметь одинаковые проценты срабатываний и тождеств при выравнивании относительно исходной последовательности. Счет от выравнивания является лучшим показателем значимости выравнивания. Этот пример использует те же гены и белки, связанные с болезнью Tay-Sachs, проанализированные в Aligning Pairs of Sequences.

Доступ к данным NCBI из рабочей области MATLAB ®

В этом примере вы будете работать непосредственно с данными о белке, так что используйте getgenpept вместо getgenbank для загрузки данных с сайта NCBI. Сначала прочитайте информацию о белке человека в MATLAB ®.

humanProtein = getgenpept('NP_000511');

Результаты поиска BLASTX, выполненного с этой последовательностью, показали, что белок Drosophila, число присоединения GenPept AAM29423, имеет некоторое сходство с последовательностью HEXA человека. Использование getgenpept чтобы загрузить эту последовательность.

flyProtein = getgenpept('AAM29423');

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

load('flyandhumanproteins.mat','humanProtein','flyProtein')
seqdisp(humanProtein)
seqdisp(flyProtein)

ans =

  10x70 char array

    '>gi|189181666|gb|NP_000511.2| beta-hexosaminidase subunit alpha pre...'
    '  1  MTSSRLWFSL LLAAAFAGRA TALWPWPQNF QTSDQRYVLY PNNFQFQYDV SSAAQPGCSV'
    ' 61  LDEAFQRYRD LLFGSGSWPR PYLTGKRHTL EKNVLVVSVV TPGCNQLPTL ESVENYTLTI'
    '121  NDDQCLLLSE TVWGALRGLE TFSQLVWKSA EGTFFINKTE IEDFPRFPHR GLLLDTSRHY'
    '181  LPLSSILDTL DVMAYNKLNV FHWHLVDDPS FPYESFTFPE LMRKGSYNPV THIYTAQDVK'
    '241  EVIEYARLRG IRVLAEFDTP GHTLSWGPGI PGLLTPCYSG SEPSGTFGPV NPSLNNTYEF'
    '301  MSTFFLEVSS VFPDFYLHLG GDEVDFTCWK SNPEIQDFMR KKGFGEDFKQ LESFYIQTLL'
    '361  DIVSSYGKGY VVWQEVFDNK VKIQPDTIIQ VWREDIPVNY MKELELVTKA GFRALLSAPW'
    '421  YLNRISYGPD WKDFYIVEPL AFEGTPEQKA LVIGGEACMW GEYVDNTNLV PRLWPRAGAV'
    '481  AERLWSNKLT SDLTFAYERL SHFRCELLRR GVQAQPLNVG FCEQEFEQT            '


ans =

  12x70 char array

    '>gi|21064387|gb|AAM29423.1| RE17456p [Drosophila melanogaster].       '
    '  1  MSLAVSLRRA LLVLLTGAIF ILTVLYWNQG VTKAQAYNEA LERPHSHHDA SGFPIPVEKS'
    ' 61  WTYKCENDRC MRVGHHGKSA KRVSFISCSM TCGDVNIWPH PTQKFLLSSQ THSFSVEDVQ'
    '121  LHVDTAHREV RKQLQLAFDW FLKDLRLIQR LDYGGSSSEP TVSESSSKSR HHADLEPAAT'
    '181  LFGATFGVKK AGDLTSVQVK ISVLKSGDLN FSLDNDETYQ LSTQTEGHRL QVEIIANSYF'
    '241  GARHGLSTLQ QLIWFDDEDH LLHTYANSKV KDAPKFRYRG LMLDTSRHFF SVESIKRTIV'
    '301  GMGLAKMNRF HWHLTDAQSF PYISRYYPEL AVHGAYSESE TYSEQDVREV AEFAKIYGVQ'
    '361  VIPEIDAPAH AGNGWDWGPK RGMGELAMCI NQQPWSFYCG EPPCGQLNPK NNYTYLILQR'
    '421  IYEELLQHTG PTDFFHLGGD EVNLDCWAQY FNDTDLRGLW CDFMLQAMAR LKLANNGVAP'
    '481  KHVAVWSSAL TNTKRLPNSQ FTVQVWGGST WQENYDLLDN GYNVIFSHVD AWYLDCGFGS'
    '541  WRATGDAACA QYRTWQNVYK HRPWERMRLD KKRKKQVLGG EVCMWTEQVD ENQLDNRLWP'
    '601  RTAALAERLW TDPSDDHDMD IVPPDVFRRI SLFRNRLVEL GIRAEALFPK YCAQNPGECI'

Первое сравнение и глобальное выравнивание

Первое, что нужно сделать - использовать seqdotplot чтобы увидеть, есть ли какие-либо области, которые четко выровнены. Это не показывает явных выравниваний, но есть некоторые области интереса.

seqdotplot(humanProtein,flyProtein,3,2)
title('Dot Plot of Two HexA-like Proteins');
ylabel('Human Protein');xlabel('Drosophila Protein');

Заметьте, что на точечном графике есть несколько диагональных отрезков. Это не особенно убедительно свидетельствует о значительном глобальном выравнивании, но можно попробовать глобальное выравнивание с помощью функции nwalign. Матрица оценки BLOSUM50 используется по умолчанию.

[sc50,globAlig50] = nwalign(humanProtein,flyProtein)
fprintf('Score = %g \n',sc50)

sc50 =

   49.6667


globAlig50 =

  3x670 char array

    'MT-S-S--R----LW----F--SLL-----LA-A-AF--A-GR------ATAL-WP----W--P-QN---FQT-----SDQR--Y---------V-LYPN---NF----Q---FQY-DVS---SAAQPGC-SVLDEAFQRY-RD--L---L-F-GSGSWPR-PYLTGK-R-HT-LE-KNVLV-VSV-V-TPG--CN-Q-----LPT--LE-SVEN---YTL-TIND-D--QCLLLSETVWGALRGLETFSQLVWKSAEGTFF--INKTEIEDFPRFPHRGLLLDTSRHYLPLSSILDTLDVMAYNKLNVFHWHLVDDPSFPYESFTFPELMRKGSYNPVTHIYTAQDVKEVIEYARLRGIRVLAEFDTPGH--T-LSWGP--GIPGLLTPCYSGSEP-S---G--TFGPVNPSLNNTYEFMSTFFLE-VSSVFP-DFYLHLGGDEVDFTCWKSNPEIQDFMRKKGFGEDFKQLESFYIQTLLDIVSSYGKGYVVWQEVFDNKVKIQPDT-I-IQVWREDI-PVNY--MKE-LELV-TKAGFRAL-LS-APWY-LNRISYGP--DWKDFYIVEPLA-FEGTPEQKALVIGGEACMWGEYVDNTNLVPRLWPRAGAVAERLWSNKLTS-DLTF----AYERLSHFRCELLRRGVQAQPLNVGFCEQE-FEQT'
    '|: : |  |    |     |  ::|     :: | |:  |  |      |:::  |    |    :|   :::     | :|  :         | ::|:   :|    |   |:  ||:   ::|:    : |: ||: : :|  |   | : ||:| |     ::| | |: ||   :|  ::  |   |   : |     | :  |: |::|   | | | ::    |  ::::: :|| :|| |::||:| : |  ::    :::::| |:| :|||:||||||:: : ||  |:  |:  |:| |||||:|  |||| |  :|||  :|:|:  :: |: |||:|| |:|:: |::|: |:|:|:|  :  :|||  |: | |: | : ::| |   |    | :||: | || ::: :: | :: : | ||: |||||||:: || :  :  | :|  |:  || :|:: : :  |   :   |  :||: :: |  :: |:: : :|||  :    ||  : :  ::: :::    |  : : |   :  : :    |:: |  :|   ::   ::|  |:|||:||| | ||:::|  |||||::|:|||||::   : |: :    :::|:| || :|:: |::|: |   :| |:  |  '
    'MSLAVSLRRALLVLLTGAIFILTVLYWNQGVTKAQAYNEALERPHSHHDASGFPIPVEKSWTYKCENDRCMRVGHHGKSAKRVSFISCSMTCGDVNIWPHPTQKFLLSSQTHSFSVEDVQLHVDTAHREVRKQLQLAFDWFLKDLRLIQRLDYGGSSSEPTVSESSSKSRHHADLEPAATLFGATFGVKKAGDLTSVQVKISVLKSGDLNFSLDNDETYQLSTQTEGHRLQVEIIANSYFGARHGLSTLQQLIWFDDEDHLLHTYANSKVKDAPKFRYRGLMLDTSRHFFSVESIKRTIVGMGLAKMNRFHWHLTDAQSFPYISRYYPELAVHGAYSE-SETYSEQDVREVAEFAKIYGVQVIPEIDAPAHAGNGWDWGPKRGM-GELAMCIN-QQPWSFYCGEPPCGQLNPKNNYTYLILQRIYEELLQHTGPTDFF-HLGGDEVNLDCWAQYFNDTD-LR--GLWCDF-MLQA-MARLKLANNGVAPKHVAVWSSALTNTKRL-PNSQFTVQVWGGSTWQENYDLLDNGYNVIFSHVDAWYLDCGFGSWRATGDAACAQYRTWQNVYKHRPWERMRLDKKRKKQVLGGEVCMWTEQVDENQLDNRLWPRTAALAERLWTDPSDDHDMDIVPPDVFRRISLFRNRLVELGIRAEALFPKYCAQNPGECI'

Score = 49.6667

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

[sc30,globAlig30] = nwalign(humanProtein,flyProtein,'scoringmatrix','blosum30')
fprintf('Score = %g \n',sc30)

sc30 =

    82


globAlig30 =

  3x670 char array

    'MT-S-S--R-----L-W--F--S-LL----L--AAAF--A-GR------ATAL-WP----W--P-QN-F--QT-----SDQR--Y---------V-LY--PN-NF----Q---F-----QY--DVS-SAAQPGCS-VLDEAFQ--RY-RDLLF-GSGSWP-RPYLTGK-R-HT-LEK-NVLV-VSV-VTP-G-CN--QLP-T-LESVE-NYTLTINDD--QC-L-L----L-SETV---W-GALRGLETFSQLVWKSAEGTFF-I-NKTEIEDFPRFPHRGLLLDTSRHYLPLSSILDTLDVMAYNKLNVFHWHLVDDPSFPYESFTFPELMRKGSYNPVTHIYTAQDVKEVIEYARLRGIRVLAEFDTPGH--T-LSWGPGIP-GLLTPCYSGSEP-S---G--TFGPVNPSLNNTYEFMSTFFLE-VSSVFP-DFYLHLGGDEVDFTCWKSNPEIQDFMRKKGFGEDFKQLESFYIQTLLDIVSSYGKGYVVWQEVFDNKVKIQPDT--IIQVWREDI-PVNY--MKE-LELV-T--KAGFRALLSAPWY-L-NRI-SYGPDWKDFYIVEPLAFEGT-PEQKALVIGGEACMWGEYVDNTNLVPRLWPRAGAVAERLWSNK-LTSDL-TF---AYERLSHFRCELLRRGVQAQPLNVGFCEQ-EFEQT'
    '|: : |  |     | :  |  : |:    :  |:|:  |  |      |:::  |    |  : :|    ::     |::|  :         | ::  |: :|    |   |     |:  |::  :::: :  ::|  ::  |: : | : ||:| |     ::| | |: ||: ::|: ::: | : | ::  |:: : |:| : |::|  ||:  |: : :    |  |::   : ||  ||:|: ||:| ::|  :: : ::::::| |:| :|||:||||||:: ::||  |:  |:::|:| |||||:|::|||| |  :|||:  |:|:: :::|::|||:||:|:|:::|::|: |:|:|:|  :  :|||    | |: |  : :| |   |  : |::||: | || ::  :: | :  : | ||: |||||||:: || : : : | :  :|:::||: |::  ::: |:: :   | ::|| ::: | :|  |::  ::|||  :: : ||  ::  : :: :  :|:: ::  : |: : : : :    |:::|   | :   : :::|: |:|||:||| | || : |  |||||::|:|||||:::    |: ::   :: |:| ||  |:  |::|: |   :|:| : | :'
    'MSLAVSLRRALLVLLTGAIFILTVLYWNQGVTKAQAYNEALERPHSHHDASGFPIPVEKSWTYKCENDRCMRVGHHGKSAKRVSFISCSMTCGDVNIWPHPTQKFLLSSQTHSFSVEDVQLHVDTAHREVRKQLQLAFDWFLKDLRLIQRLDYGGSSSEPTVSESSSKSRHHADLEPAATLFGATFGVKKAGDLTSVQVKISVLKSGDLNFSLD-NDETYQLSTQTEGHRLQVEIIANSYFGARHGLSTLQQLIWFDDEDHLLHTYANSKVKDAPKFRYRGLMLDTSRHFFSVESIKRTIVGMGLAKMNRFHWHLTDAQSFPYISRYYPELAVHGAYSE-SETYSEQDVREVAEFAKIYGVQVIPEIDAPAHAGNGWDWGPKRGMGELAMC-INQQPWSFYCGEPPCGQLNPKNNYTYLILQRIYEELLQHTGPTDFF-HLGGDEVNLDCW-A-QYFND-TDLRGLWCDFM-LQA-MARLKLANNGVAPKHVAVWSSALTN-TKRLPNSQFTVQVWGGSTWQENYDLLDNGYNVIFSHVDAWYLDCGFGSWRATGDAACAQYRTWQNVYKHRPWERMRLDKKRKKQVLGGEVCMWTEQVDENQLDNRLWPRTAALAERLWTDPSDDHDMDIVPPDVFRRISLFRNRLVELGIRAEALFPKYCAQNPGECI'

Score = 82

Это дает выравнивание, которое имеет некоторые области довольно сильного сходства, но является ли это выравнивание статистически значимым? Один из способов выяснить, является ли этот счет значительным - использовать методы Монте-Карло. Учитывая, что последовательность мухи была найдена с помощью поиска BLAST, существуют некоторые доказательства того, что существует сходство между этими двумя последовательностями. Разумно ожидать, что счет для такого выравнивания будет выше, чем счета, полученная при выравнивании случайных последовательностей аминокислот к белку.

Оценка значимости счета

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

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

rng(0,'twister')
n = 50;
globalscores = zeros(n,1);
flyLen = length(flyProtein.Sequence);
for i = 1:n
    perm = randperm(flyLen);
    permutedSequence = flyProtein.Sequence(perm);
    globalscores(i) = nwalign(humanProtein,permutedSequence,'scoringmatrix','blosum30');
end

Теперь постройте график счетов как штриховой график. Обратите внимание, что поскольку вы используете случайным образом сгенерированные последовательности.

figure
buckets = ceil(n/5);
hist(globalscores,buckets)
hold on;
stem(sc30,1,'k')
title('Determining Alignment Significance using Monte Carlo Techniques');
xlabel('Score'); ylabel('Number of Sequences');

Счета выравниваний к случайным последовательностям могут быть аппроксимированы распределением экстремальных значений типа 1. Используйте evfit функция из Statistics and Machine Learning Toolbox™, чтобы оценить параметры этого распределения.

parmhat = evfit(globalscores)

parmhat =

  -31.7597    6.6440

Наложите график функции плотности вероятностей предполагаемого распределения.

x = min(globalscores):max([globalscores;sc30]);
y = evpdf(x,parmhat(1),parmhat(2));
[v, c] = hist(globalscores,buckets);
binWidth = c(2) - c(1);
scaleFactor = n*binWidth;
plot(x,scaleFactor*y,'r');
hold off;

Из этого графика можно увидеть, что глобальное выравнивание (globAlig30) явно статистически значимо.

Пример, где счет не статистически значима

На веб-сайте FLYBASE вы можете найти все гены дрозофилы бета-N-ацетилгексосаминидазы. Ген, на который вы смотрели до сих пор, упоминается как CG8824. Теперь вы хотите взглянуть на другой подобный ген, например Hexo1.

flyHexo1 = getgenpept('AAL28566');

Последовательность fly Hexo1 аминокислот также предусмотрена в MAT-файле flyandhumanproteins.mat.

load('flyandhumanproteins.mat','flyHexo1')
seqdisp(humanProtein)

ans =

  10x70 char array

    '>gi|189181666|gb|NP_000511.2| beta-hexosaminidase subunit alpha pre...'
    '  1  MTSSRLWFSL LLAAAFAGRA TALWPWPQNF QTSDQRYVLY PNNFQFQYDV SSAAQPGCSV'
    ' 61  LDEAFQRYRD LLFGSGSWPR PYLTGKRHTL EKNVLVVSVV TPGCNQLPTL ESVENYTLTI'
    '121  NDDQCLLLSE TVWGALRGLE TFSQLVWKSA EGTFFINKTE IEDFPRFPHR GLLLDTSRHY'
    '181  LPLSSILDTL DVMAYNKLNV FHWHLVDDPS FPYESFTFPE LMRKGSYNPV THIYTAQDVK'
    '241  EVIEYARLRG IRVLAEFDTP GHTLSWGPGI PGLLTPCYSG SEPSGTFGPV NPSLNNTYEF'
    '301  MSTFFLEVSS VFPDFYLHLG GDEVDFTCWK SNPEIQDFMR KKGFGEDFKQ LESFYIQTLL'
    '361  DIVSSYGKGY VVWQEVFDNK VKIQPDTIIQ VWREDIPVNY MKELELVTKA GFRALLSAPW'
    '421  YLNRISYGPD WKDFYIVEPL AFEGTPEQKA LVIGGEACMW GEYVDNTNLV PRLWPRAGAV'
    '481  AERLWSNKLT SDLTFAYERL SHFRCELLRR GVQAQPLNVG FCEQEFEQT            '

Повторите процесс генерации глобального выравнивания и затем используя случайные сочетания аминокислот, чтобы оценить значимость глобального выравнивания.

[Hexo1score,Hexo1Alignment] = nwalign(humanProtein,flyHexo1,'scoringmatrix','blosum30')
fprintf('Score = %g \n',Hexo1score)
Hexo1globalscores = zeros(n,1);
flyLen = length(flyHexo1.Sequence);
for i = 1:n
    perm = randperm(flyLen);
    permutedSequence = flyHexo1.Sequence(perm);
    Hexo1globalscores(i) = nwalign(humanProtein,permutedSequence,'scoringmatrix','blosum30');
end

Hexo1score =

  -72.2000


Hexo1Alignment =

  3x534 char array

    'MTSSRL-WFSLLLAAAFA-GRATALWPWPQNFQTSDQRYVLYPNNFQFQYDVSSAAQPGCSVLDEAFQRYRDLLFGSGSWPRPYLTGKRHTLEKNVLVVSVVTPGC-NQLPTLESVENYTLTINDDQCLLLSETVWGALRGLETFSQLVWKSAEGTFFINKTEIEDFPRFPHRGLLLDTSRHYLPLSSILDTLDVMAYNKLNVFHWHLVDDPSFPYESFTFPELMRKGSYNPVTHIYTAQDVKEVIEYARLRGIRVLAEFDTPGHTLSWGPGIPGLLTPCYSGSEPSGTFGPVNPSLNNTYEFMSTFFLEVSSVFPDFYLHLGGD-EV-DFTCWKSNPEIQDFMRKKGFGEDFKQLESFYIQTLLDIVSSYGKGYVVWQEVFDNKVKIQPDTIIQVWREDIPVNYMKELELVTKAGFRALLSAPWYLNRISYGPDWKDFYIVEPLAFEGTPEQKALVIGGEACMWGEYVDNTNLVPRLWPRAGAVAERLWSNKLTSDLTFAYERLSHFRCELLRRGVQAQPLNVGFCEQEFEQT'
    '|:  :|  |   :::: :   ::: :|  ::: : :|| | |  : : : ||::::: |  |   :: |::   | : : | ::  | :|   || :     | :| |: |  :|   : : : : :|  |::|| :::  :::: : ::    || | ::  ::|:  | |:|   |   : :  :|  :: :  :  :  : :| : ::::: :|:  |   || ::: : |::: :: |::|:  : :  | :|  : | |  |:  | |  :|:|::   :  :| | |:|::::::| :: : : ||: : ::|| :||:  | | : | | | |  : :|  :  | : |:|  |:  ::    :|:  ::| | :| :  :  |:  ::|    |  :: : :| ::   |  |   |  | :: |  |:    :|   :  :: |: || :|   : |  :: : : |: |    : :|   : ::  |   |        |:   ||     |:| :: : :: :'
    'MALVKLNTFHWHITDSHSFPLEVKKRPELHKLGAYSQRQV-Y--T-R-R-DVAEVVEYG-RV--RGI-RVMP-EF-D-A-PAHVGEGWQH---KN-M-----T-ACFNAQP-WKS---F-C-V-EPPCGQLDPTV-NEM--YDVL-EDIY----GTMF-DQF-NPDI--F-HMG--GD---E-VS-TSCWNS-S--Q--P--IQQW-M-KKQGWGLETADF---MRLWGHFQ-TEAL-GR-VDKVANGTHT-PI-IL--W-TSG--LTEEPFIDEYLNPERYIIQ-IWTTG-VDPKVKKILE-RG-YKIIVSN-YDALYLDCGGAGWVTDGNNWCS-PYI-GW-QKV-Y--D-NSLKS--IAGDYEH-HVLGAEGAIWSEQID-EHTL--DN--RFW----P--RA-S-AL-AE---R--L---W-SNPAE-G--WR--Q-AES-RLL-LHRQR-LVDNG---L-G--AE-A-MQPQ-W-CL-Q-NE-H-ECPI--D---A--------CS---RGSGRLGLIVLLLLTTLS-A'

Score = -72.2

Постройте графики счетов, вычислите параметры распределения и наложите PDF на столбчатую диаграмму.

figure
buckets = ceil(n/5);
hist(Hexo1globalscores,buckets)
title('Determining Alignment Significance using Monte Carlo Techniques');
xlabel('Score');
ylabel('Number of Sequences');
hold on;
stem(Hexo1score,1,'c')
parmhat = evfit(Hexo1globalscores)
x = min(Hexo1globalscores):max([Hexo1globalscores;Hexo1score]);
y = evpdf(x,parmhat(1),parmhat(2));
[v, c] = hist(Hexo1globalscores,buckets);
binWidth = c(2) - c(1);
scaleFactor = n*binWidth;
plot(x,scaleFactor*y,'r');
hold off;

parmhat =

  -70.6926    7.0619

В этом случае кажется, что выравнивание не является статистически значимым. Более высокие выравнивания оценки могут быть легко сгенерированы из случайного сочетания аминокислот в последовательности. Можно вычислить приблизительное p-значение из предполагаемого экстремального значения CDF: Однако для получения достоверной оценки параметров крайнего значения pdf, из которых можно вычислить достаточно точное p-значение, необходимо намного больше 50 случайных сочетание.

p = 1 - evcdf(Hexo1score,parmhat(1),parmhat(2))

p =

    0.4458

Одна вещь, которая нужно заметить, что длины двух последовательностей очень различаются. Длина HEXA1 человека составляет 529 остатков, а длина белка Hexo1 мухи составляет всего 383 остатка. Когда вы пытаетесь выровнять эти две последовательности глобально, это различие в длине означает, что в последовательность должно быть введено большое количество погрешностей. Это означает, что значимость счетов будет сильно зависеть от GAPOPEN и EXTENDGP параметры. (См. справку для nwalign для получения дополнительной информации Вместо использования глобального выравнивания в этом случае лучшим подходом может быть просмотр локального выравнивания между этими двумя последовательностями.

Использование локальной трассы и Randseq

Теперь вы повторите процесс оценки значимости выравнивания на этот раз с помощью локального выравнивания и немного другого метода генерации случайных последовательностей. Вместо простой перестановки букв в последовательности, альтернативой является получение последовательности из полиномиального распределения, которое оценивается из последовательности белка мухи. Вы можете сделать это с помощью aacount и randseq функции; первый оценивает аминокислотные частоты последовательности запроса, а последующий случайным образом создает новые последовательности на основе этого распределения.

[lscore,locAlig] = swalign(humanProtein,flyHexo1,'scoringmatrix','blosum30')
fprintf('Score = %g \n',lscore)

localscores = zeros(n,1);
aas = aacount(flyHexo1);
for i = 1:n
    randProtein = randseq(flyLen,'FROMSTRUCTURE',aas);
    localscores(i) = swalign(humanProtein,randProtein,'scoringmatrix','blosum30');
end

lscore =

   152


locAlig =

  3x361 char array

    'MAYNKLNVFHWHLVDDPSFPYESFTFPELMRKGSYNPVTHIYTAQDVKEVIEYARLRGIRVLAEFDTPGHT-LSWG-PGIPGLL-TPCYSGSEPSGTFGPVNPSLNNTYEFMSTFFLEVSSVF-PDFYLHLGGDEVDFTCWKSNPEIQDFMRKKGFG-E--DFKQLES-FYIQTL--LD-IVSSYGKGYVVWQE-VFDNK-V-K-IQPD-TIIQVWREDI-P-VNYMKELEL-VTKAGFRALLSAPWYLNRISYGPDWKDFYI-VEP-L--AFEG-TPEQKALVIGGEACMWGEYVDNTNLVPRLWPRAGAVAERLWSNKLTSDLTFAYER-LSHFRCELLRRGVQAQPLNVGFCEQ-EFE'
    '||: |||:||||::|::|||:|    |||:: |:|::  ::|| :||:||:||:|:|||||: |||:|:|:  :|  ::::::: :: ::: : : : |:::|::| :|:::: ::  : : | ||:: |:|||||: :||:|:::|| :|:|:|:| |  ||::| : | :::|  :| ::::     ::|:: : : : : : : |:  |||:|   : | |: : |  : :: ::: ||: : :  :::: | :| : ||  :: :  :::: : : :  |:|:|: :|:| :|  :|  |:||||:|:|||||||: ::    | :| | | |  |:  |: |: :   :| | | |'
    'MALVKLNTFHWHITDSHSFPLEVKKRPELHKLGAYSQR-QVYTRRDVAEVVEYGRVRGIRVMPEFDAPAHVGEGWQHKNMTACFNAQPWKSFCVEPPCGQLDPTVNEMYDVLEDIYGTMFDQFNPDIF-HMGGDEVSTSCWNSSQPIQQWMKKQGWGLETADFMRLWGHFQTEALGRVDKVANGTHTPIILWTSGLTEEPFIDEYLNPERYIIQIWTTGVDPKVKKILERGYKIIVSNYDALYLDCGGAGWVTDGNNWCSPYIGWQKVYDNSLKSIAGDYEHHVLGAEGAIWSEQIDEHTLDNRFWPRASALAERLWSNP-AEGWRQAESRLLLH-RQRLVDNGLGAEAMQPQWCLQNEHE'

Score = 152

Постройте графики счетов, вычислите параметры распределения и наложите PDF на столбчатую диаграмму.

figure
hist(localscores,buckets)
title('Determining Alignment Significance using Monte Carlo Techniques');
xlabel('Score');
ylabel('Number of Sequences');
hold on;
stem(lscore,1,'r')
parmhat = evfit(localscores)
x = min(localscores):max([localscores;lscore]);
y = evpdf(x,parmhat(1),parmhat(2));
[v, c] = hist(localscores,buckets);
binWidth = c(2) - c(1);
scaleFactor = n*binWidth;
plot(x,scaleFactor*y,'r');
hold off;

parmhat =

   40.8331    3.9312

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

При локальном выравнивании кажется, что выравнивание статистически значимо. Фактически, просмотр локального выравнивания показывает очень хорошее выравнивание для всей длины Hexo1 последовательности.

close all;

Документация по Bioinformatics Toolbox

Поддержка