Uniprot整合数据比对分析

1.从uniprot上下载蛋白序列以及带有物种分类的excel文件

首先对excel文件进行处理:将文件变为以下形式

image.png
然后对fasta文件进行处理
下载的蛋白序列如下：
image.png
在比对后映射的时候只需要“|”中间的蛋白号，所以用python对序列进行预处理。

import sys

with open(sys.argv[1]) as f:
        with open(sys.argv[2],"a") as g:
                for i in f.readlines():
                        if(i.startswith(">")):
                                i = i.strip("\n").split("|")
                                g.write(">"+str(i[1])+"\n")
                        else:
                                g.write(str(i))

处理完数据如下：

image.png

2.利用diamond比对

diamond makedb --in end_p_gene.fasta --db end_p_gene
diamond blastx --db end_p_gene.dmnd -q unigene.fa -o end_p_dna.m8 -t /dev/shm &

3.对diamond结果进行处理：筛选p值小于0.00001的结果；然后根据预测基因ID去冗余。

awk -F"\t" '$11<0.00001' p_dna.m8 > p_dna_1.m8
sort -k1,1 -u p_dna_1.m8 > p_dna_2.m8

4.用R语言对结果进行整理(以下例举两个基因)

映射基因和物种

p_gene <- read.table("p_dna_2.m8",sep = "\t")
appa <- read.table("appa.csv",sep = ",",header = T)
gcd <- read.table("gcd.csv",sep = ",",header = T)
mergi <- function(a,b,c){
  colnames(a)[2] <- c
  f <- merge(a,b,by = c)
}
appa_1 <- mergi(a= p_gene , b= appa , c= "Entry")
gcd_1 <- mergi(a= p_gene , b= gcd , c= "Entry")

结果如下：

image.png

整合丰度表并将各基因数据合并

ab <- read.table("unigene_cov.xls",sep = "\t",header = T)
name <- list(appa_1,gcd_1)
p_gene_ab <- data.frame()
for(i in 1:2){
    colnames(name[[i]])[2] <- "Gene_ID"
    name[[i]] <- name[[i]][,-(3:12)]
    name[[i]] <- merge(name[[i]],ab,by = "Gene_ID")
    name[[i]] <- name[[i]][,c(1:10,(ncol(name[[i]])-35):ncol(name[[i]]))]
    p_gene_ab <- rbind(name[[i]],p_gene_ab)
}
write.csv(p_gene_ab,"p_gene_ab.csv",row.names = F)

结果如下：

image.png
image.png

若想对单个基因进行物种贡献分析
以下例举Family水平上的物种对基因丰度的贡献度
注意：跑完千万要把x和endgame格式化

appa_1ab <- name[[1]]
gcd_1ab <- name[[2]]
x <- x[,c(6:8,(ncol(x)-35):ncol(x))]
x <- x[order(x$Family),]
endgame <- data.frame()
for(i in 2:nrow(x)){
  if(i == nrow(x) & x$Family[i] == x$Family[i-1]){
    for (e in (ncol(x)-35):ncol(x)) {
      x[i,e] <- sum(x[i-1,e],x[i,e])
    }
    endgame <- rbind(as.data.frame(x[i,]),endgame)
  }else if(i == nrow(x) & x$Family[i] != x$Family[i-1]){
    endgame <- rbind(as.data.frame(x[i,]),endgame)
  }else if(i != nrow(x) & x$Family[i] == x$Family[i-1]){
    for (e in (ncol(x)-35):ncol(x)){
      x[i,e] <- sum(x[i-1,e],x[i,e])
    }
  }else{
    endgame <- rbind(as.data.frame(x[i-1,]),endgame)
  } 
}
gcd_Family_2 <- endgame

结果如下：

image.png
注释：对单个基因的物种贡献分析代码也可以对每种基因的丰度进行累加，得到基因丰度。