本文主要是介绍NGS项目五：用R语言生成极端嗜盐古菌蛋白序列进化树，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！
基因表达数据分析和实验设计密不可分，总体来说，实验设计有两大类思路：一类是时间序列分析，主要思想是测定基因多个时间点的表达值，通过聚类和主成分分析等分析手段寻找调控基因，进而研究其深层机制；第二类是基因表达差异的限制性分析。
所谓GO注释，就是将表示基因或其产物的ID映射到一组GO的ID上，用这组GOterm来描述这个基因。而实际应用中，人们更关心一组基因（比如差异表达基因）的共同点，因此会对它们所对应的所有GO的分布情况进行分析，有利于发现新的现象，或者集中到感兴趣的方向。
GO富集分析的统计基础是超几何分布，简单来说就是根据下列Fisher精确检验（Fisherexact test)公式,对每个GOterm计算一个P值：
P=(M/k)((N-M)/(n-k))/(N/n)
通路（Pathway)分析包括通路注释和通路富集分析。通路富集分析的基本思路、统计模型等基本和GO富集分析一致。常用的公共通路数据库主要有KEGG、BioCarta和GenMAPP。由于KEGG注释来源和授权问题，Biocondocutor已经转而使用更加开源的Reactome数据库。
基因组或转录组序列的注释分析包括基因识别和基因功能注释两个方面。基因识别的核心是确定全基因组序列中所有基因的确切位置。主要有两种方法：一是通过相似性比对一直基因和蛋白序列得到简接证据；二是基于各种统计模型和算法从头预测。对从头预测出的基因进行高通量功能注释也是基于相似性搜索。
ID映射可以通过编程自动化进行，大大提高了数据注释的效率。除了自行编程，常用的数据库ID之间的映射可以使用Uniprot数据库网站的工具。芯片数据处理过程中，经常将Affymetrix芯片组的ID映射到其他数据库ID,这方面的常用工具是NIH网站的DAVID(http://david.abcc.ncifcrf.gov/)。
对于一组差异基因，研究人员若关心这组基因富集到哪些基因本体论术语或者通路上，这就需要应用统计学方法来对富集的显著性进行定量。统计为后续研究指明了方向，是生物信息学研究的核心；而可视化可以帮组研究人员更好地理解统计结果，提供进一步研究的思路和灵感。R语言提供了强大的一体化的统计和绘图功能。
反复阅读、调试源代码是学习计算机编程语言的唯一捷径。
极端嗜盐古菌发展出许多复杂且完整的分子机制以对抗盐分、氧气、光和养分的变动，因此可以尝试从蛋白水平寻找某些序列特征，以期找到一些相关的机制。可分为以下五个方面：
A.将蛋白质序列批量导入;
B.统计每条序列的氨基酸百分比含量;
C特定模序匹配与统计，提取包含模序2次以上的蛋白质序列;
D 对提取的序列进行两两比对;
E 根据比对结果构建系统发育树。
A.将蛋白质序列批量导入;
>seq_import<-function(input_file){
+my_fasta<-readLines(input_file);
+y<-regexpr("^>",my_fasta,perl=T);
+ y[y==1]<-0;
+index<-which(y==0);
+distance<-data.frame(start=index[1:(length(index)-1)],end=index[2:length(index)]);
+distance<-rbind(distance,c(distance[length(distance[,1]),2],length(y)+1));
+distance<-data.frame(distance,dist=distance[,2]-distance[,1]);
+seq_no<-1:length(y[y==0]);
+index<-rep(seq_no,as.vector(distance[,3]));
+my_fasta<-data.frame(index,y,my_fasta);
+my_fasta[my_fasta[,2]==0,1]<-0;
+seqs<-tapply(as.vector(my_fasta[,3]),factor(my_fasta[,1]),paste,collapse="",simplify=F);
+seqs<-as.character(seqs[2:length(seqs)]);
+Desc<-as.vector(my_fasta[c(grep("^>",as.character(my_fasta[,3]),perl=TRUE)),3]);
+my_fasta<-data.frame(Desc,Length=nchar(seqs),seqs);
+Acc<-gsub(".*gb\\|(.*)\\|.*","\\1",as.character(my_fasta[,1],perl=T);
+my_fasta<-data.frame(Acc,my_fasta);
+my_fasta;
}
B.统计每条序列的氨基酸百分比含量;
pattern_match<-function(pattern,sequences,hit_num){
pos<-gregexpr(pattern,as,character(sequences[,4]),perl=T);
posv<-unlist(lapply(pos,paste,collapse=“，”））;
posv[posv==-1]<-0;
hitsv<-unlist(lapply(pos,function(x)if(x[1]==-1){0}else{length(x)}));
sequences<-data.frame(sequences[,1:3],Position=as.vector(posv),His=hitsv,sequences[,4]);
tag<-gsub(“([A-Z]”,“\\L\\1”,as.character(sequences[sequences[,5]>hit_num,6]),perl=T,ignore.case=T);
pattern2=paste(“(“,pattern,”)”,seq=””);
tag<-gsub(pattern2,“\\U\\,1”,tag,perl=T,ignore.case=T);
export<-data.frame(sequences[sequences[sequences[,5]>hit_num,-6],tag);
export<-data.frame(Acc=paste(“>”,export[,1],seq=“”),seq=export[,6];
write.table(as.vector(as.character(t(export))),file=“Hit_sequences.fasta”,quote=F,row.names=F,col.names=F);
cat(“含有模序\””,pattern,”\”超过”，hit_num,“个的所有蛋白质序列已经写入当前工作目录下文件“Hit_sequence.fasta””,”\n”,seq=””);
selected<-sequences[sequences[,5]>hit_num,];
cat(“极端嗜盐古菌蛋白组中以下序列含有模序\””,pattern,”\”的数量超过2个：”，“\n”,seq=””);
print(selected[,1:5]);
selected;
}
C特定模序匹配与统计，提取包含模序2次以上的蛋白质序列;
getAApercentage<-function(sequences){
AA<-data.frame(AA=c(“A”,”C”,”D”,”E”,”F”,”G”,”H”,”I”,”K”,”L”,”M”,”N”,”P”,”Q”,”R”,”S”,”T”,”V”,”W”,”Y”));
Aastat<-lapply(strsplit(as.character(sequences[,6]),””),table);
for(i in1:length(AAstat)){
Aaperc<-Aastat[[i]]/sequences[,3][i]*100;
Aaperc<-as.data.frame(AAperc);
names(AAperc)[2]<-as.vector(sequences[i,1]);
AA<-merge(AA,Aaperc,by.x=”AA”,by.y=”Varl”,all=T);
}
for(i in1:length(AA[[1]])){
for（jin 1:length(AA)){
if(is.na(AA[i,j])){
AA[i,j]<-0;
}
}
}
Aapercentage<-data.frame(AA,Mean=apply(AA[,2:length(AA)],1,mean,na.rm=T));
write.csv(AApercentage,file=”AApercentage.csv”,row.names=F,quote=F);
AApercentage;
}
D 对提取的序列进行两两比对;
seq_alignment<-function(sequences){
shell(“del/fmy_needle_file”);
for(i in1:length(sequences[,1])){
cat(as.character(paste(“>”,as.vector(sequences[i,1]),seq=””)),as.character(as.vector(sequences[i,6])),file=”file1”,seq=”\n”);
for(j in1:length(sequences[,1])){
cat(as.character(paste(“>”,as.vector(sequences[j,1]),seq=””)),as.character(as.vector(sequences[i,6])),file=”file1”,seq=”\n”);
shell(“needlefile1 file2 stdout -gapopen 10.0 -gapextend 0.5 >>my_needle_file\”\n”);
}
}
cat(“Needle程序完成所有序列的两两比对，结果存入文件\”my_needle_file\”\n”);
}
E1定义函数’getScoreMatrix‘求得得分矩阵
getScoreMatrix<-function(sequences){
score<-readLine(“my_needle_file”);
score<-score[grep(“^#Score”,score,perl=T)];
score<-gsub(“.*”,””,as.character(score),perl=T);
score<-as.numeric(score);
scorem<-matrix(score,length(sequences[,1]),length(sequence[,1]),dimnames=list(as.vector(sequences[,1]),as.vector(sequence[,1])));
scorem.dist<-as.dist(1/scorem);
hc<-hclust(scorem.dist,method=”complete”);
plot(hc,hang=-1,main=”DistanceTree Based on Needle All-Against-All Comparison”,xlab=”sequencename”,ylab=”distance”);
scorem;
}
E2定义函数’infile_produce'
infile_produce<-function(scorem){
z<-l/scorem;
len=sqrt(length(scorem));
z[seq(1,length(scorem),by=(len+1))]<-0;
z<-round(z,7);
write.table(len,file=”infile”,quote=F,row.names=F,col.names=F);
write.table(as.data.frame(z),file=”infile”,append=T,quote=F,col.names=F,seq=”\t”);
cat(“Phylip格式的距离矩阵已经输出到工作目录下名为‘infile'的文件，以便使用Phylip软件继续进行分析。”，“\n”);
}
课题实现
A调用函数’seq_import'导入数据
setwd(“C:/workingdirectory”);
my_file<-”AE004437.faa”;
my_sequence<-seq_import(input_file=my_file);
B调用函数‘pattern_match'寻找模序
hit_sequences<-pattern_match(pattern=”H..H{1,2}”,sequences=my_sequences,hit_num=2)
C调用函数’getAApercentage'统计氨基酸百分含量
AA_percentage<-getAApercentage(sequences=hit_sequences)
D调用函数‘seq_alignment'进行序列两两比对
seq_alignment(sequences=hit_sequences)
E1调用函数’getScoreMatrix'得到得分矩阵
score_matrix<-getScoreMatrix(sequences=hit_sequences)
E2调用函数‘infile_produce'生成PHYLIP软件的输入文件infile
infile_produce(scorem=score_matrix)
E3 调用PHYLIP软件生产进化树
将infile文件拖入neighbor窗口
E4 进化树查看与编辑
使用retree.exe
这篇关于NGS项目五：用R语言生成极端嗜盐古菌蛋白序列进化树的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！
