华侨大学 - 生物信息学 - 实验报告 2

华侨大学计算机科学与技术学院

生物信息学实验报告

姓名：_______纪新婷_________ 学号：_1025112018_____ 专业班级：_计算机科学与技术1班________

1.matlab软件Bioinformatics Toolbox工具包中的常用基本函数： （1）baselookup('Complement', SeqNT)

baselookup('Complement', 'ｕｇｇｇａｕｕａｕｃｃａａａａｕｇｕｇａ') ans =

ｕｇｇｇａｕｕａｕｃｃａａａａｕｇｕｇａ ans =

ｕｇｇｇａｕｕａｕｃｃａａａａｕｇｕｇａ

函数 baselookup()的第一个叁数表示的是该函数所使用的方法，在这个例子

中的“Complement”表示的是该函数所实现的功能是对核苷酸序列的补充，而”SeqNT”表示的是核苷酸的序列，在这个例子中我们引用得非典的核苷酸序列。

（2）SeqRNA = dna2rna(SeqDNA)

SeqRNA = dna2rna('CCTGCCTTAATGTGCCTC')

SeqRNA =

CCUGCCUUAAUGUGCCUC

函数dna2rna(SeqDNA)实现的功能是将DNA序列转化为RNA序列，函数中的参数SeqDNA是DNA序列。

（3）SeqR = seqreverse(SeqNT) s = 'CCTGCCTTAATGTGCCTC';

seqreverse(s)

ans =

CTCCGTGTAATTCCGTCC

函数seqreverse(s)返回一个DNA的核苷酸序列的反向链，函数中的参数s是DNA序

列。

（4）gag = multialignread('aagag.aln'); showalignment(gag)

函数showalignment(gag)的作用是显示多序列比对。

（5）ramachandran(PDBid) ramachandran('1E31')

函数ramachandran(PDBid)的功能是蛋白质数据库（PDB）的数据绘制的图， 其中‘1E13’是人的Survivin基因在PDB中的标识符。

（6）seqdotplot(Seq1, Seq2)

seqdotplot('CCTGCCTTAATGTGCCTC','CTCCGTGTAATTCCGTCC')

函数seqdotplot(Seq1, Seq2)的作用是创建两个序列的散点图，以达到可视化的两个序列之间的匹配。

2.学习Bioinformatics Toolbox工具包中的seqstatsdemo.m和aligndemo.m两个例子：

(1)seqstatsdemo.m:

% 大家公认的人类线粒体基因组序列的编号为：NC_012920，我们仅仅使用核苷酸

（nucleotide）序列，所以用'SequenceOnly'这个标志来限制。利用getgenbank函数就能将我们需要的人类线粒体基因组的核苷酸序列读入到matlab的工作空间。

mitochondria = getgenbank('NC_012920','SequenceOnly',true); % 如果你没有连接网络，您可以从MAT-file使用commandload线粒体，取消这个如果没有互联网连接

%matlab中的‘whos’的信息时给出序列的大小 whos mitochondria

%nt前缀表示核苷酸，使用ntdensity 函数可以得到该核苷酸的组成 ntdensity(mitochondria) 其结果如下：

% ntdensity 函数得到的图已经表明该核苷酸富含A-T，但是我们能使用basecount函数

得到更具体的信息：A、T、C、G的数目

ans =

A: 5113 C: 5192 G: 2180 T: 4086

>> basecount(seqrcomplement(mitochondria)) % basecount函数得到的是5'--3'单链的A、T、C、G的数目，而 seqrcomplement函数能得到其补链（3'-5'）。

ans =

A: 4086 C: 2180 G: 5192 T: 5113

>> basecount(mitochondria,'chart','pie');

% 将A、T、C、G的数目以饼状图的形式表示出来。

>> dimercount(mitochondria,'chart','bar')

% 将核苷酸序列中的二聚体（dimer）的信息即数目显示出来。二聚体：相邻两个核苷酸的组合。

ans =

AA: 1594 AC: 1495 AG: 801 AT: 1223 CA: 1536 CC: 1779 CG: 439 CT: 1438 GA: 615 GC: 716 GG: 427 GT: 421 TA: 1368 TC: 1202 TG: 512 TT: 1004

>> codoncount(mitochondria)

% 将核苷酸中的密码子（codon）的数目显示出来，codonconut函数没有使用选项，则用第一种阅读框架显示出来（the first reading frame）。密码子：由三个相邻的核苷酸组成的信使核糖核酸（mRNA）基本编码单位，共有64种密码子，有61钟氨基酸密（amino acid）码子和3种终止密码子。在蛋白质合成时，一种氨基酸密码子代表一种氨基酸。

AAA - 172 AAC - 157 AAG - 67 AAT - 123 ACA - 153 ACC - 163 ACG - 42 ACT - 130 AGA - 58 AGC - 90 AGG - 50 AGT - 43 ATA - 132 ATC - 103 ATG - 57 ATT - 96 CAA - 166 CAC - 167 CAG - 68 CAT - 135 CCA - 146 CCC - 215 CCG - 50 CCT - 182 CGA - 33 CGC - 60 CGG - 18 CGT - 20 CTA - 187 CTC - 126 CTG - 52 CTT - 98 GAA - 68 GAC - 62 GAG - 47 GAT - 39 GCA - 67 GCC - 87 GCG - 23 GCT - 61

GGA - 53 GGC - 61 GGG - 23 GGT - 25 GTA - 61 GTC - 49 GTG - 26 GTT - 36 TAA - 136 TAC - 127 TAG - 82 TAT - 107 TCA - 143 TCC - 126 TCG - 37 TCT - 103 TGA - 64 TGC - 35 TGG - 27 TGT - 25 TTA - 115 TTC - 113 TTG - 37 TTT - 99 >> for frame = 1:3 figure subplot(2,1,1); codons{1,frame} = codoncount(mitochondria,'frame',frame,'figure',true); title(sprintf('Codons for Frame %d of Human Mitochondria Genome',frame)); subplot(2,1,2); codons{2,frame} = codoncount(mitochondria,'reverse',true,'frame',frame,'figure',true); title(sprintf('Codons for Reverse Frame %d of Human Mitochondria Genome',frame));

End

% 使用codoncount函数的其他的3种阅读框架显示出来。

Exploring the Open Reading Frames (ORFs)

%在核苷酸序列中明显可以看出，如果有任何的开放阅读框架。一个ORF序列的DNA或RNA，可以把bepotentially翻译成蛋白质。* seqshoworfs *功能是可用于可视化的ORF序列中。

%注意：在HTML教程只显示输出的第一页，然而当运行演示，你将能够检查完成项目基因组中使用的滚动条图。

seqshoworfs(mitochondria);

%如果你比较这个输出显示在NCBI网页似乎有较少的ORFs基因，比预期的基因还少。脊椎动物线粒体不使用标准遗传密码，一些密码子的线粒体基因组中不同的含义。为了在MATLAB中能更多的了解不同的遗传密码，那就查看功能*遗传密码的帮助。 help geneticcode

>> orfs = seqshoworfs(mitochondria,'GeneticCode','Vertebrate Mitochondrial',... 'alternativestart',true)

% 以脊椎动物的基因编码方式显示核苷酸序列的ORFs，注意到在第一种阅读框架下，这多出两个ORFs,一个起始位置在4471，另一个在5905.他们刚好对应 ND2 (NADH

dehydrogenase subunit 2) 和COX1 (cytochrome c oxidase subunit I) 基因。

orfs =

1x3 struct array with fields: Start Stop

Extracting and analyzing the ND2 protein

ORF起始位置在4471的基因对应的蛋白质就是ND2 protein. ND2Start = 4471;

startIndex = find(orfs(1).Start == ND2Start) ND2Stop = orfs(1).Stop(startIndex) 找到ND2蛋白质的终止密码子。

startIndex =

24

ND2Stop =

5512

ND2Seq = mitochondria(ND2Start:ND2Stop);

知道ND2的终止密码子之后就能得到ND2对应的编码区域 codoncount(ND2Seq)

利用前面的codoncount函数就能显示出该序列的密码子。

>>

%

>> % >> %

>> aminolookup('code',nt2aa('CTA'))

aminolookup('code',nt2aa('ATC'))

% 如果不记得密码子对应的氨基酸，则可以用nt2aa和aminolookup函数来查找。

ans =

Leu Leucine

ans =

Ile Isoleucine

>> ND2 = nt2aa(ND2Seq,'GeneticCode','Vertebrate Mitochondria');

% nt2aa函数将核苷酸序列转换成对应的氨基酸序列。值得注意的是，这儿的基因编码方式必须限定在vertebrate mitochondrial genetic code。 figure

>> aacount(ND2,'chart','bar')

title('Histogram of Amino Acid Count for Human Mitochondrial Genome');

% 使用account函数得到完整的氨基酸组成。

%通知高亮氨酸，苏氨酸和异亮氨酸含量也缺乏半胱氨酸或天冬氨酸。

%%你可以使用* * * atomiccomp molweight功能去发现更多的关于ND2蛋白。 ND2AtomicComp = atomiccomp(ND2); ND2MolWeight = molweight(ND2);

>> proteinplot(ND2)

% 使用proteinplot函数得到一个分析蛋白质性质的GUI，能够轻松得到蛋白质的性质，比如：疏水性（hydrophobicity）。

%%你也可以创建使用* proteinpropplot *序列各性能曲线。

proteinpropplot(ND2,'PropertyTitle','Parallel beta strand')

Inspecting Other Annotated Features

%%你也可以看看所有的已注释的人类线粒体基因组的特点。下载完整的GenBank条目”nc_001807和* featuresparse *功能，探讨其特点。

mitochondria_gbk = getgenbank('NC_001807'); features = featuresparse(mitochondria_gbk); genes = {features.gene.gene};

ND2annotations = features.gene(10); % ND2 is in the 10th position

features =???? source: [1x1 struct]?? D_loop: [1x1 struct]?? gene: [1x37 struct]?? tRNA: [1x22 struct]?? rRNA: [1x2 struct]?? STS: [1x28 struct]?? misc_feature: [1x1 struct]?? CDS: [1x13 struct]??????coding_sequences_id =????ND1 ND2 COX1 COX2 ATP8 ATP6 COX3 ND3 ND4L ND4 ND5 ND6 CYTB????

%%创建一个地图显示所有使用* featuresmap *功能在GenBank进入发现的特点。 [h,l] = featuresmap(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')

%%* <地址：bioinfo-feedback@mathworks.com？个体反馈% 20 % % % % 20seqstatsdemo 20

20bioinformatics 20toolbox % 202.5表明该演示增强> *。 displayEndOfDemoMessage(mfilename)

（2）Aligndemo.m： Aligning pairs of sequences

在基因库中提取出核苷酸序列，找到其ORFs后，需要使用全局和局部的比对算法来比对（Align）这些序列。

Using the Help Browser as a web browser

>> web('http://www.ncbi.nlm.nih.gov/') % 访问NCBI网站。

>> web('http://www.ncbi.nlm.nih.gov/books/bv.fcgi?call=bv.View..ShowSection&rid=gnd') % NCBI中的\部分非常有价值，这部分对医学遗传学（medical genetics）作了一个全面的介绍。 >> web('http://www.ncbi.nlm.nih.gov/books/bv.fcgi?call=bv.View..ShowSection&rid=gnd.section.220') >> web('http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?val=13128865&db=Nucleotide&dopt=GenBank')

% 这一部分我们将以Tay-Sachs Disease的基因为例。Tay-Sachs是一种由15号染色体上两种等位的基因（HEXA）突变引起的常染色体隐性遗传病。NCBI中HEXA的编号是：NM_000520。

%?cessing NCBI Data from the MATLAB(R) Workspace

%你可以使用MATLAB中的* getgenbank *函数读取序列信息。

humanHEXA = getgenbank('NM_000520');

%%通过BLAST搜索或搜索在小鼠基因组你会发现ak080777是鼠标己糖胺酶基因

mouseHEXA = getgenbank('AK080777');

%%如果你没有一个活的网络连接，您可以从一个MAT文件加载数据使用命令负荷< = =取消导致如果没有连接

Exploring the Open Reading Frames (ORFs)

%%你可以使用这个seqshoworfs * *功能寻找人类HEXA基因序列中的ORFs。注意，最长ORF在第一个阅读框。在变量humanorfs的输出值中在每个阅读框架的所有基因密码子是一个结构的启动位置和停止。 humanORFs = seqshoworfs(humanHEXA.Sequence); humanORFs =

1x3 struct array with fields: Start

Stop

%%看看现在的老鼠六基因的ORF。在这种情况下，ORF为框架1。 mouseORFs = seqshoworfs(mouseHEXA.Sequence); mouseORFs =

1x3 struct array with fields: Start

Stop

Aligning the Sequences

%%第一步是使用全局序列比对，寻找这些序列之间的相似性。你可以看看核苷

酸序列之间的对准，但它通常是更有益的看蛋白序列之间的对准。所以，你要做的第一件事是将核苷酸序列转换成相应的氨基酸序列。使用* nt2aa *函数来做这个。

humanProtein = nt2aa(humanHEXA.Sequence); mouseProtein = nt2aa(mouseHEXA.Sequence);

%%一个寻找序列之间的相似性的最简单的方法是用一个点阵图。

seqdotplot(humanProtein,mouseProtein) ylabel('Human hexosaminidase A'); xlabel('Mouse hexosaminidase A');

%%用默认的设置，点图有点难以解释的，所以你可以稍微更严格的点阵图。

seqdotplot(humanProtein,mouseProtein,4,3)

ylabel('Human hexosaminidase A');xlabel('Mouse hexosaminidase A');

%%对角线表明有可能是一个很好的定位，你可以看看现在的全局对准使用功能* nwalign *使用Needleman Wunsch算法。

[~, globalAlignment] = nwalign(humanProtein,mouseProtein);

%在帮助浏览器* *功能* showalignment *匹配和相似的残基在不同的颜色突出显示对齐。

showalignment(globalAlignment);

Refining the Alignment

%对准很好除了终端段。请注意，对齐的区域划定的开始（m-methionine）和停止（*）序列中的氨基酸。如果序列缩短，只有翻译区域被认为是，那么可能你会得到一个更好的对准。使用*找到命令来寻找在每个序列的起始氨基酸指数：

humanStart = find(humanProtein == 'M',1); mouseStart = find(mouseProtein == 'M',1);

%%同样，使用*到*命令查找索引的第一站的翻译开始后发生的。需要特别注意因为还停留在* humanprotein *序列的开始。

humanStop = find(humanProtein(humanStart:end)=='*',1) + humanStart - 1;

mouseStop = find(mouseProtein(mouseStart:end)=='*',1) + mouseStart - 1;

humanStart = 70

mouseStart = 11

humanStop = 599

mouseStop = 539

%%利用这些指标对截形序列。

humanSeq = humanProtein(humanStart:humanStop); humanSeqFormatted = seqdisp(humanSeq);

mouseSeq = mouseProtein(mouseStart:mouseStop); mouseSeqFormatted = seqdisp(mouseSeq); humanSeqFormatted =

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*

mouseSeqFormatted =

1 MAGCRLWVSL LLAAALACLA TALWPWPQYI QTYHRRYTLY PNNFQFRYHV SSAAQAGCVV

61 LDEAFRRYRN LLFGSGSWPR PSFSNKQQTL GKNILVVSVV TAECNEFPNL

ESVENYTLTI

121 NDDQCLLASE TVWGALRGLE TFSQLVWKSA EGTFFINKTK IKDFPRFPHR GVLLDTSRHY

181 LPLSSILDTL DVMAYNKFNV FHWHLVDDSS FPYESFTFPE LTRKGSFNPV THIYTAQDVK

241 EVIEYARLRG IRVLAEFDTP GHTLSWGPGA PGLLTPCYSG SHLSGTFGPV NPSLNSTYDF

301 MSTLFLEISS VFPDFYLHLG GDEVDFTCWK SNPNIQAFMK KKGFTDFKQL ESFYIQTLLD

361 IVSDYDKGYV VWQEVFDNKV KVRPDTIIQV WREEMPVEYM LEMQDITRAG FRALLSAPWY

421 LNRVKYGPDW KDMYKVEPLA FHGTPEQKAL VIGGEACMWG EYVDSTNLVP RLWPRAGAVA

481 ERLWSSNLTT NIDFAFKRLS HFRCELVRRG IQAQPISVGC CEQEFEQT*

%%如果你看这两个序列，你就会看到一个很好的全局对准。

[~, alignment] = nwalign(humanSeq,mouseSeq); showalignment(alignment);

%%或者，您可以使用编码区的信息（CDS）从GenBank中的数据结构来发现基因的编码区。你也可以从这个结构编码区中翻译。 humanCodingRegion = humanHEXA.Sequence(...

humanHEXA.CDS.indices(1):humanHEXA.CDS.indices(2)); mouseCodingRegion = mouseHEXA.Sequence(...

mouseHEXA.CDS.indices(1):mouseHEXA.CDS.indices(2))； humanTranslatedRegion = humanHEXA.CDS.translation; mouseTranslatedRegion = mouseHEXA.CDS.translation;

Local Alignment

%%而不是截断序列寻找更好的调整，一个 替代的方法是使用一个局部比对。功能* SWALIGN *执行使用Smith Waterman算法局部比对。这表明整个编码区和几个残基在两末端的基因相似的一个很好的定位。

[~, localAlignment] = swalign(humanProtein,mouseProtein); showalignment(localAlignment);

Alignment of Complementary DNA Sequences

%%多序列比对的补充而不是身份认同是一个用户定义的得分矩阵。在这种情况下，1分是给予补充和1则。评分矩阵，定义以下的行和列的顺序为：C G碱基T.从鼠标六基因第一30个核苷酸，将与它的补充。

scoring_matrix = [-1,-1,-1,1;-1,-1,1,-1;-1,1,-1,-1;1,-1,-1,-1]; [score, compAlignment] = nwalign(mouseHEXA.Sequence(1:30), ... seqcomplement(mouseHEXA.Sequence(1:30)), 'SCORINGMATRIX', ... scoring_matrix, 'ALPHABET', 'NT'); scoring_matrix =

-1 -1 -1 1 -1 -1 1 -1 -1 1 -1 -1 1 -1 -1 -1 score = 30

compAlignment =

GCTGCTGGAAGGGGAGCTGGCCGGTGGGCC :::::::::::::::::::::::::::::: CGACGACCTTCCCCTCGACCGGCCACCCGG

本文来源：https://www.bwwdw.com/article/047p.html