diamond比对以及结果数据处理

利用diamond比对数据

注意：diamond在无root权限的linux运行比对命令时需要加参数-t /dev/shm创建临时文件

#建库--in nr输入文件，--db nr 输出的数据库前缀
diamond makedb --in nr --db nr
#比对
diamond blastx --db nr -q reads.fna -o dna_fmt6
diamond blastp --db nr -q reads.faa -o protein_fmt6

注意事项:
默认参数主要是针对段短序列，对于比较长的序列，使用--sensitive或--more-senstive提高敏感度。
默认的e-value阈值是0.001, 而BLAST是10，因此会比BLAST结果更加严格
性能优化:
设置比较低的-e参数
设置-k参数，减少输出的联配数目。这会降低临时文件大小和最终结果
--top会输出得分比最好的分数低一定百分比的结果，
--compress 1: 输出结果会以gzip进行压缩
————————————————
原文链接：https://blog.csdn.net/u012110870/article/details/102804629
————————————————
有时间可以查查blastx和blastp的原理，它们会把预测的基因按照不同的数据框（应该是6种）读取然后比对
所以比对结果有许多重复的预测基因id，对应相同的数据库基因id和不同的E值，分数等。

（diamond结果和blast结果排列格式一样，有12列。以NCycDB为例）

因此我们需要对结果进行处理：

1. 根据E值筛选（这个也可以在diamond参数里调整）一般是取小于0.0001或0.00001（1e-5）

#第11列是E值,这里根据0.00001筛选
awk -F"\t" '$11<0.00001' dna_fmt6 > dna_fmt6_1e-5

2. 筛选分数最高得预测基因id的结果

sort -k1,1 -u dna_fmt6_1e-5 > dna_fmt6_1e-5_one

以下处理工具为R语言

3. 将结果的id和id_基因名的映射文件（map文件）结合。

e <- read.table(" dna_fmt6_1e-5_one",sep = "\t")
f <- read.table("id2gene.map",sep = "\t")
colnames(e)[2] <- "pi"
colnames(f)[1] <- "pi"
g<-merge(e,f,by="pi")

4. 将结果结合丰度表

u <- read.table("unigene_cov.xls",sep = "\t",header = T)
result<-merge(g,u,by="Gene_ID")
#result1<- plyr::join(g,u,by="Gene_ID") 

注意：当用merge（）函数时会将缺失值NA忽略。而join（）函数会将有缺失值NA的基因也匹配。
所以我觉得用merge（）就好了

5. 将相同基因名（也就是同源基因）的丰度相加得到最终结果。

最新的优化代码在《物种丰度处理》那一章，快很多倍。

x <- result
x <-x[order(x$Gene_ID),]
ending <- as.data.frame(t(as.data.frame(c(colnames(x)))))
colnames(ending) <- colnames(x)
for(i in 2:length(x$Gene_ID)){
  if(i == length(x$Gene_ID) & x$Gene_ID[i] == x$Gene_ID[i-1]){
    for(e in 2:37) {
        x[i,e] <- sum(x[i-1,e],x[i,e])
    }
    ending <- rbind(as.data.frame(x[i,]),ending)
  }else if(i == length(x$Gene_ID) & x$Gene_ID[i] != x$Gene_ID[i-1]){
    ending <- rbind(as.data.frame(x[i,]),ending)
  }else if(i != length(x$Gene_ID) & x$Gene_ID[i] == x$Gene_ID[i-1]){
    for (e in 2:37){
      x[i,e] <- sum(x[i-1,e],x[i,e])
    }
  }else{
    ending <- rbind(as.data.frame(x[i-1,]),ending)
  } 
}
write.table(ending,file="ending.csv",sep=',',quote = F,row.names=F)