评价真菌引物覆盖度(5)R计算引物覆盖度

在得到Blastn的结果后,即可结合数据库的物种信息在R中计算引物的覆盖度。注意Blastn的输出格式为7。以UNITE数据库,计算phylum水平上的覆盖度为例:

统计unite数据库phylum水平上的物种信息

>database = read.table(file="total_taxa_unite.txt",sep="\t")

#统计总序列数量

>seqnum = nrow(database)

#统计phylum的个数。第6列是phylum的信息。

>dphy=unique(database[,6])

Blastn结果中匹配上的phylum

这里计算Mismatch=0时的覆盖度。其条件为mismatch=0,identify=100%,aligh length = primer length。

(x[,13]=="100")&

(x[,14]==primerlength)&

(x[,15]==mismatch))

#第7列为phylum信息。统计匹配上的门的个数

>phylum=unique(x_mis0[,7])

计算覆盖度

###1.计算匹配的序列占数据库序列的百分比。第三列为序列id,用于统计匹配上的序列数。

#注意,这个地方要取unique。否则有简并存在序列会重复计算。覆盖度可能会大于100%。

>matchseq = length(unique(x_mis0[,3]))

>total_coverage = paste("total coverage is",matchseq*100/seqnum,"%")

###2.phylum水平上覆盖度

#匹配上的门

>intphy = intersect(phylum,dphy)##交集

>intphy =c("intersectphylum number is", length(intphy),intphy)

>phylum_coverage = paste("phylum coverage is",length(intphy)*100/length(dphy),"%")

#没有匹配上的门

>difphy = setdiff(dphy,phylum)##差集

>difphy =c("diffphylum number is", length(difphy),difphy)

最后输出total_coverage 和 phylum_coverage,intphy,difphy

>sink("coverage_phylum.txt")

>total_coverage

>phylum_coverage

>intphy

>difphy

>sink()

系列历史

一个环境工程专业却做生信分析的深井冰博士,深受拖延症的困扰。想给自己一点压力,争取能够不定期分享学到的生信小技能,亦或看文献过程中的一些笔记与小收获,记录生活中的杂七杂八。

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20181206G0X1OJ00?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 yunjia_community@tencent.com 删除。

扫码关注云+社区

领取腾讯云代金券