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

Доступ к данным NCBI из MATLAB® Workspace

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

humanProtein = getgenpept('NP_000511');

Результаты поиска BLASTX, выполняемого с этой последовательностью, показали, что белок Дрозофилы, инвентарный номер 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)
showalignment(globAlig50);

Score = 49.6667 

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

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

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 можно искать всю Дрозофилу beta-N-acetylhexosaminidase гены. На ген, на который вы смотрели до сих пор, ссылаются как 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)
showalignment(Hexo1Alignment);
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

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: Однако намного больше, чем 50 случайных перестановок необходимы, чтобы получить надежную оценку экстремума параметры PDF, от которых можно вычислить довольно точное p-значение.

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

p =

    0.4458

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

Используя Local Alignment и Randseq

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

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

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

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.