代码比较复杂~~~~还没理解透彻
library(VennDiagram)
library(plotrix)
library(RColorBrewer)
library(scales)
data<-read.table('taxa.txt',row.names = 1,header=T,comment.char ='',sep='\t')
#导入OTU数据
meta<-read.table('group.txt',row.names = 1,header=T,comment.char ='',sep='\t')
#导入分组数据
taxonomy<-data[,dim(data)[2]]
#提取最后一列到taxonomy
data<-data[,-dim(data)[2]]
#去掉最后一列
meta<-meta[match(colnames(data),rownames(meta)),]
#将分组和OTU数据根据样品名match
data<-data[,!is.na(meta)]
meta<-meta[!is.na(meta)]
abdc<-colSums(t(data))
func1<-function(x){
return(tapply(x,INDEX = meta,sum)>0)
}
tb1<-apply(data,1,func1)
unig<-row.names(tb1)
dl<-colSums(tb1)
ll<-length(unig)
com<-sum(dl==4)
com_spe<-data.frame(OTUID=colnames(tb1)[dl==4],taxonomy=taxonomy[dl==4],sum_abundance=abdc[dl==4])
write.table(com_spe,"common_species.xls",sep = "\t",row.names = F)
for (i in 1:4) {
is_uni<-(dl==1&tb1[i,]==1)
flower_data[i]<-sum(is_uni)
uni_spe<-data.frame(OTUID=colnames(tb1)[is_uni],taxonomy=taxonomy[is_uni],sum_abundance=abdc[is_uni])
write.table(uni_spe,paste(out,"/",group,"_",unig[i],"_special_species.xls",row.names = F,sep="\t")
}
is_uni<-(dl==1&tb1[4,]==1)
flower_data[i]<-sum(is_uni)
uni_spe<-data.frame(OTUID=colnames(tb1)[is_uni],taxonomy=taxonomy[is_uni],sum_abundance=abdc[is_uni])
write.table(uni_spe,"4_special_species.xls",row.names = F,sep="\t")
if(length(unig)<=5)
ls1<-list()
for (i in 1:4) {ls1[[i]]<-colnames(tb1)[tb1[i,]]}
names(ls1)<-rownames(tb1)
venn.diagram(ls1,filename = 'output.tiff',
imagetype="tiff",
alpha=0.50,
lwd=1.2,
cat.cex=1.4,
col = 'black',
fill=rainbow(length(ls1)),margin=0.15)
venn.diagram(
ls1,
filename = 'output3.tiff',
col = 'black',
fill = brewer_pal(palette = 'Set3')(4)
)
