annot.R 3.3 KB
Newer Older
1
2
#!/usr/bin/env RScript

3
args=commandArgs(trailingOnly=TRUE)
4
5
#ddsRdata=args[1]
OUTDIR=args[1]
6
7
cond1=args[2]
cond2=args[3]
8
9
ASSEMBLY<-args[4]
CORRANNOT<-args[5]
10
11
12

library(clusterProfiler)

13
14
comp<-paste(cond1,cond2,sep="__vs__")

15
16
17
18
19
20
makeGeneList<-function(data,org){
	names=data[,1]
	s=strsplit(as.character(names),"[|]")
	symb=do.call(rbind,s)
	data[,1]=symb[,2]

21
22
23
24
	pAdj=data[,"padj"]
	InvPAdj=1-pAdj
	data=data.frame(data, InvPAdj)
	logC=data[,"log2FoldChange"]
25
26
	logC[logC<0] <- -1
	logC[logC>=0] <- 1
27
	data$InvPAdj=data$InvPAdj*logC
28
29
30
31
	geneID=bitr(data$Name, fromType="SYMBOL",toType=c("ENTREZID"),OrgDb=org)
	doublons=which(duplicated(geneID$SYMBOL))
	if(length(doublons)>0){geneID=geneID[-doublons,]}
	data=data[is.element(data$Name, geneID$SYMBOL),]
32
	geneList=data[,dim(data)[2]]
33
34
35
36
37
38
	names(geneList)=as.character(geneID$ENTREZID)
	geneList=sort(geneList, decreasing=TRUE)		
	return(geneList)
}

# Récupération des identifiants de l'organisme pour GO et KEGG (annotations fonctionnelles) dans le fichier de correspondance assemblage->annotations
Eric CHARPENTIER's avatar
Eric CHARPENTIER committed
39
corresIDorg=read.table(CORRANNOT,header=T,row.names=1,stringsAsFactors=T)
40
orgGO=as.character(corresIDorg[ASSEMBLY,1])
41
42
orgKegg=as.character(corresIDorg[ASSEMBLY,2])

Eric CHARPENTIER's avatar
Eric CHARPENTIER committed
43
44
45
46
47
#if((!require(orgGO,character.only=TRUE))){
#	requireNamespace("BiocManager", quietly = TRUE)
#	BiocManager::install(orgGO, updates = TRUE)
#	require(orgGO,character.only=TRUE)
#}
48

Eric CHARPENTIER's avatar
Eric CHARPENTIER committed
49
data=read.csv(paste(paste(OUTDIR,comp,sep="/"),"DEseqResFiltered.tsv",sep="/"), header=TRUE, sep="\t",stringsAsFactors=T)
50

51
52
#Sys.setenv(http_proxy="http://cache.ha.univ-nantes.fr:3128")
#Sys.setenv(https_proxy="https://cache.ha.univ-nantes.fr:3128")
53

54
55
56
57
if(dim(data)[1]>10){
	geneList=makeGeneList(data,orgGO)
	gene=names(geneList)
	ego=enrichGO(gene=gene,OrgDb=orgGO,ont="ALL",pAdjustMethod = "BH",pvalueCutoff  = 0.01,qvalueCutoff  = 0.05,readable      = TRUE)
58
	png(paste(paste(OUTDIR,comp,sep="/"),"dotplotGO.png",sep="/"),width=1000,height=600)
59
60
61
	if(dim(ego)[1] != 0){
		print(clusterProfiler::dotplot(ego))
	}
62
	dev.off()
63
	write.table(ego,paste(paste(OUTDIR,comp,sep="/"),"annotGo.tsv",sep="/"),sep="\t",row.names=F,quote=F)
64
65

	ekegg=enrichKEGG(gene, organism = as.character(orgKegg), keyType = "kegg", pvalueCutoff = 0.05,pAdjustMethod = "BH", minGSSize = 10, maxGSSize = 500,qvalueCutoff = 0.2, use_internal_data = FALSE)
66
	eKeggSymbol=setReadable(ekegg,orgGO,keyType="ENTREZID")
67
	png(paste(paste(OUTDIR,comp,sep="/"),"dotplotKEGG.png",sep="/"),width=1000,height=600)
68
69
70
	if(dim(eKeggSymbol)[1] != 0){
		print(clusterProfiler::dotplot(eKeggSymbol))
	}
71
	dev.off()
72
	write.table(eKeggSymbol,paste(paste(OUTDIR,comp,sep="/"),"annotKegg.tsv",sep="/"),sep="\t",row.names=F,quote=F)
73
74
}

Eric CHARPENTIER's avatar
Eric CHARPENTIER committed
75
dataTot=read.csv(paste(paste(OUTDIR,comp,sep="/"),"DEseqRes.tsv",sep="/"), header=TRUE, sep="\t",stringsAsFactors=T)
76
77
geneListTot=makeGeneList(dataTot,orgGO)

Eric CHARPENTIER's avatar
Eric CHARPENTIER committed
78
ego2=gseGO(geneList=geneListTot,OrgDb=get(orgGO),ont="ALL",nPerm=1000,minGSSize=20,maxGSSize=500,pvalueCutoff=1,verbose=FALSE)
79
ego2Symbol=setReadable(ego2,orgGO,keyType="ENTREZID")
80
81

kk2=gseKEGG(geneList=geneListTot,organism=as.character(orgKegg),nPerm=1000,minGSSize=2,pvalueCutoff=1,verbose=FALSE)
82
kk2Symbol=setReadable(kk2,orgGO,keyType="ENTREZID")
83

84
85
write.table(ego2Symbol,paste(paste(OUTDIR,comp,sep="/"),"gseGo.txt",sep="/"),sep="\t",row.names=F,quote=F)
write.table(kk2Symbol,paste(paste(OUTDIR,comp,sep="/"),"gseKegg.txt",sep="/"),sep="\t",row.names=F,quote=F)