我用维斯塔对两种细菌进行了基因组比较。
这个工具给出了两个细菌之间常见的DNA序列区域,但我最感兴趣的是知道在一个缺乏第二个的细菌中存在哪些CDS。
通过使用R,我设法使用VISTA信息生成一个data.frame,其中包括第一个细菌所独有的碱基区域(范围)。这些区域必须含有第二个缺失的基因(CDS's)。
head(rango_vacio) # Regions (mapped bp) exclusive to the first bacteria
V1 V2
11552 13259
13365 13263
37168 37169
..... .....另一方面,我处理了同样细菌的gff文件来提取CDS序列。此数据包含每个CDS的开始和结束,以及相应蛋白质的登录名。
head(cds_TIGR4) # A list of the cds of this bacteria
startbp endbp accession
197 1158 NP_344444
1717 2853 NP_344445
2864 3112 NP_344446
..... .... .....重要的是:数据帧"rango_vacio“和"cds_TIGR4”都使用与基本相同的引用,所以我可以比较两者。
现在,我的问题的答案应该很容易完成,因为我只需要使用CDS本身的范围来找出每个rango_vacio范围中存在哪些CDS
我可以通过使用一组非常复杂的for循环来完成这一任务,但是我想向你们中的任何人学习这是否可以通过任何其他更短的方法来完成。
发布于 2015-03-26 22:55:11
最后,我相信我已经找到了自己的方法
GenomicRanges不能在我的情况下使用,因为我的data.frame中有一个包含cds范围,包括strandness。另一个只包含一个范围
所以我使用了IRange包。
我简化了两个数据帧,以包含范围和cds的开始和结束。一个叫兰戈,另一个叫cds
library(IRanges)
ir_rango <- IRanges(rango[,1], rango[,2])
ir_cds <- IRanges(cds[,1], cds[,2])
common <- findOverlaps(ir_cds, ir_rango)
common <- as.matrix(common)
unique_cds <- common[,1]
uniques <- which(duplicated(unique_cds))
uniquesuniques包含ir_cds中显示的相应范围的行号。现在我只需要提取cds的名称
https://stackoverflow.com/questions/29278028
复制相似问题