В этом примере показан метод, который можно использовать для исследования значимости выравнивания последовательностей. Количество удостоверений или положительных результатов в трассе не является четким показателем значительной трассы. Перестановка последовательности из выравнивания будет иметь аналогичные проценты положительных результатов и идентичности при выравнивании относительно исходной последовательности. Оценка от выравнивания является лучшим показателем значимости выравнивания. В этом примере используются те же гены и белки, связанные с болезнью Тэя-Сакса, проанализированные в выравнивающих парах последовательностей.
В этом примере вы будете работать непосредственно с белковыми данными, чтобы использовать 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 функция из 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');
Последовательность мухи 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 для получения дополнительной информации Вместо использования глобального выравнивания, в этом случае лучшим подходом может быть просмотр локального выравнивания между двумя последовательностями.
Теперь вы повторите процесс оценки значимости выравнивания на этот раз с использованием локального выравнивания и несколько другого метода генерации случайных последовательностей. Вместо простого перестановки букв в последовательности альтернативой является получение последовательности из полиномиального распределения, которое оценивается по последовательности белка мухи. Это можно сделать с помощью 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;