文档章节

ORF识别代码

g
 gwc19880401
发布于 2015/04/02 09:28
字数 731
阅读 12
收藏 1

library(Biostrings)
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\\|(.*)\\|.*","\\l",as.character(my_fasta[,1]),perl=T)
my_fasta<-data.frame(Acc,my_fasta)
my_fasta
}
data=seq_import("AvsB NCBI seq_name.fasta")
mydata=readLines("AvsB NCBI seq_name.fasta")
seq_id=regexpr("^>",mydata,perl=TRUE)
gene_info=mydata[seq_id==1]
gi= sub("^>gi\\|([0-9]*)\\|.*","\\1",gene_info,perl=TRUE)
gene_name=sub("^>gi\\|.*\\|(.*)","\\1",gene_info,perl=TRUE)
data_vector=matrix(c(gi,gene_name,rep("",length(gi)*14)),nrow=length(gi),ncol=16)
colnames(data_vector)=c("Entrez ID","Gene_infor","ORF1","ORF2","ORF3","TC_ORF1","TC_ORF2","TC_ORF3","AA_ORF1","AA_ORF2","AA_ORF3","AA_TC_ORF1","AA_TC_ORF2","AA_TC_ORF3","ORF","Length")
for (i in 1:length(gi)){

data_vector[i,3]=as.character(data[i,4])
data_vector[i,4]=as.character(DNAString(as.character(data_vector[i,3]))[2:length(DNAString(as.character(data_vector[i,3])))])
data_vector[i,5]=as.character(DNAString(as.character(data_vector[i,3]))[3:length(DNAString(as.character(data_vector[i,3])))])

data_vector[i,6]=as.character(reverseComplement(DNAString(as.character(data_vector[i,3]))))
data_vector[i,7]=as.character(reverseComplement(DNAString(as.character(data_vector[i,4]))))
data_vector[i,8]=as.character(reverseComplement(DNAString(as.character(data_vector[i,5]))))

data_vector[i,9]=as.character(translate(DNAString(as.vector(data_vector[i,3]))))
data_vector[i,10]=as.character(translate(DNAString(as.vector(data_vector[i,4]))))
data_vector[i,11]=as.character(translate(DNAString(as.vector(data_vector[i,5]))))
data_vector[i,12]=as.character(translate(DNAString(as.vector(data_vector[i,6]))))
data_vector[i,13]=as.character(translate(DNAString(as.vector(data_vector[i,7]))))
data_vector[i,14]=as.character(translate(DNAString(as.vector(data_vector[i,8]))))
for (h in 9:14){
a=unlist(strsplit(as.character(data_vector[i,h]),"",useBytes=TRUE))=="M"

if (length(unique(a))==1){ data_vector[i,h]="FMA*K"} else {

a=unlist(strsplit(as.character(data_vector[i,h]),"",useBytes=TRUE))=="*"

if (length(unique(a))==1){

ab=as.numeric(unlist(gregexpr("M",as.character(data_vector[i,h]),useBytes=TRUE)))

data_vector[i,h]=substr(as.character(data_vector[i,h]),min(ab),nchar(as.character(data_vector[i,h])))

} else {

ab=as.numeric(unlist(gregexpr("M",as.character(data_vector[i,h]),useBytes=TRUE)))
ac=as.numeric(unlist(gregexpr("\\*",as.character(data_vector[i,h]),useBytes=TRUE)))

ah=matrix(as.numeric(""),nrow=length(ab),ncol=3)
colnames(ah)=c("AALength","l_value","k_value")

for (k in 1:length(ab)) {

ag=matrix(as.numeric(""),nrow=length(ac),ncol=3)
colnames(ag)=c("AALength","l_value","k_value")

for (l in 1:length(ac)) {

ag[l,1]=ac[l]-ab[k]
ag[l,2]=l
ag[l,3]=k  }

if (max(ag[,1])<0) {
ah[k,1]=0
ah[k,2]=1
ah[k,3]=1} else {

aj=data.frame(ag)
aj=aj[aj[,1]>0,]
aj=matrix(unlist(aj),ncol=3)
ak=min(aj[,1])

if (nrow(aj)==0) {
ah[k,1]=0
ah[k,2]=aj[m,2]
ah[k,3]=aj[m,3]

} else{

for (m in 1:nrow(aj)) {

if ( aj[m,1]==ak) {
ah[k,1]=aj[m,1]
ah[k,2]=aj[m,2]
ah[k,3]=aj[m,3]} else {}}} } }

ao=max(ah[,1])

for (n in 1:length(ab) ) {

if ( ah[n,1]==ao ) {

ap=ah[n,2]
aq=ah[n,3] } else {} }

ax=ab[aq]
ay=ac[ap]-1
data_vector[i,h]=as.character( paste(unlist(strsplit(as.character(data_vector[i,h]),"",useBytes=TRUE))[ax:ay],sep="",collapse=""))

}}}

for (j in 9:14) {
if (as.numeric(nchar(data_vector[i,j]))==max(as.numeric(nchar(data_vector[i,9])),as.numeric(nchar(data_vector[i,10])),as.numeric(nchar(data_vector[i,11])),as.numeric(nchar(data_vector[i,12])),as.numeric(nchar(data_vector[i,13])),as.numeric(nchar(data_vector[i,14]))))
{ data_vector[i,15]=as.character(data_vector[i,j])} else {}
}
}
data_vector[,16]=nchar(as.character(data_vector[,15]))
data_vector=cbind(gene_info,data_vector)
AAseq=matrix("",nrow=nrow(data_vector)*2,ncol=1)
for ( xy in 1:nrow(data_vector)) {
AAseq[2*xy-1,1]=as.character(data_vector[xy,1])
AAseq[2*xy,1]=as.character(data_vector[xy,16])
}
write.table(AAseq,"AvsB_AAseq_2.txt",sep="",row.names=FALSE,col.names=FALSE,quote=FALSE)

 

© 著作权归作者所有

g
粉丝 0
博文 1
码字总数 731
作品 0
广州
私信 提问
【译】R包介绍:Online Random Forest

作者:顾全,浙江大学软件工程硕士,现任桃树科技算法工程师 地址:https://github.com/ZJUguquan/OnlineRandomForest 参与:Cynthia 翻译:本文为天善智能编译,未经容许,禁止转载 介绍 On...

天善智能
2018/05/17
0
0
ROSALIND Bioinformatics(16-20)

刷题ROSALIND,练编程水平 http://rosalind.info/problems/list-view/ 16. Finding a Protein Motif (寻找蛋白质模体) image.png 解题思路 为表示蛋白质模体的多种形式,设计以下表示形式:[...

thinkando
2017/11/16
0
0
蓝梦图片碎片重组复原处理软件----蓝梦软件BestRecoveryForImage

蓝梦图片碎片重组复原处理软件----蓝梦软件BestRecoveryForImage 软件截图: 全面支持Canon(CR2,CRW),Nikon(NEF),Olympus(ORF),Sony(ARW,SR2),Fujifilm(RAF),Panasonic(RAW), Pentax(PEF),Si......

aijuanwok
2014/05/05
0
0
Aspose.BarCode已修复关于PDF417条码识别和生成的各种问题条码控件网

Aspose.BarCode是由Aspose Pty Ltd所开发的一款功能强大,且稳健的条形码生成和条码识别组件,其使用托管的C#编写,能帮助开发者快速简便的向其Microsoft应用程序(WinForms, ASP.NET 和.NE...

yidongkaifa
2014/08/18
180
0
OpenCV3.3人脸识别模块的API的变化

1. 前言 开始用最新版的OpenCV进行人脸识别的小伙伴也许已经发现了,人脸识别的最新API变化了。这也正是人脸识别等contrib模块没有放到主仓库的原因:不稳定,仍在开发中。当然这次的变化也预...

冰不语
2017/11/21
0
0

没有更多内容

加载失败,请刷新页面

加载更多

java通过ServerSocket与Socket实现通信

首先说一下ServerSocket与Socket. 1.ServerSocket ServerSocket是用来监听客户端Socket连接的类,如果没有连接会一直处于等待状态. ServetSocket有三个构造方法: (1) ServerSocket(int port);...

Blueeeeeee
今天
6
0
用 Sphinx 搭建博客时,如何自定义插件?

之前有不少同学看过我的个人博客(http://python-online.cn),也根据我写的教程完成了自己个人站点的搭建。 点此:使用 Python 30分钟 教你快速搭建一个博客 为防有的同学不清楚 Sphinx ,这...

王炳明
昨天
5
0
黑客之道-40本书籍助你快速入门黑客技术免费下载

场景 黑客是一个中文词语,皆源自英文hacker,随着灰鸽子的出现,灰鸽子成为了很多假借黑客名义控制他人电脑的黑客技术,于是出现了“骇客”与"黑客"分家。2012年电影频道节目中心出品的电影...

badaoliumang
昨天
15
0
很遗憾,没有一篇文章能讲清楚线程的生命周期!

(手机横屏看源码更方便) 注:java源码分析部分如无特殊说明均基于 java8 版本。 简介 大家都知道线程是有生命周期,但是彤哥可以认真负责地告诉你网上几乎没有一篇文章讲得是完全正确的。 ...

彤哥读源码
昨天
15
0
jquery--DOM操作基础

本文转载于:专业的前端网站➭jquery--DOM操作基础 元素的访问 元素属性操作 获取:attr(name);$("#my").attr("src"); 设置:attr(name,value);$("#myImg").attr("src","images/1.jpg"); ......

前端老手
昨天
7
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部