评价真菌引物覆盖度（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()

系列历史

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