annotate rmsk_summary_table_multiple.r @ 5:378565f5a875 draft

Uploaded
author petr-novak
date Fri, 22 Nov 2019 07:56:48 -0500
parents a4cd8608ef6b
children 62fefa284036
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
1 #! /usr/bin/env Rscript
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
2 # repeat masker summary
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
3 # analysis of *.out file
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
4 # input arguments: <rmsk_result.out> <reads.fas.
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
5 # calculates totoal proportion of matched repetetive families and averagee score and print the table
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
6 suppressPackageStartupMessages(library(getopt))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
7
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
8 # INPUT ARGUMENTS:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
9 opt = getopt(matrix(c(
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
10 'fasta', 'f', 1, "character",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
11 'repeat_masker','r', 1, "character",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
12 'output', 'o', 1, "character",
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
13 'help','h',0,'logical'),ncol=4,byrow=T))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
14
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
15 if (!is.null(opt$help)) {
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
16 cat("Usage:\n\n -f fasta file \n -r repeat masker output\n -o output files\n")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
17 stop()
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
18 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
19
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
20 RMfiles=system(paste("ls ",opt$repeat_masker,sep=""),intern=T)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
21 fasta=system(paste("ls ",opt$fasta,sep=""),intern=T)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
22
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
23
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
24
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
25 getsubdir=function(x){
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
26 xx=gsub(".*/","",x)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
27 sdir=gsub(xx,"",x,fixed=T)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
28 sdir
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
29 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
30
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
31 fasta.summary=function(reads){
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
32 suppressPackageStartupMessages(library(Biostrings,quiet=T))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
33 if (exists("readDNAStringSet")){
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
34 seqs=readDNAStringSet(reads)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
35 }else{
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
36 seqs=read.DNAStringSet(reads)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
37 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
38
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
39 N=length(seqs)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
40 Ls=nchar(seqs)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
41 L=sum(Ls)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
42 M=median(Ls)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
43 A=mean(Ls)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
44 output=c(N,L,M,A)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
45 names(output)=c("N","length","median","mean")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
46 output
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
47 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
48
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
49 repCont=function(repname){ # uses external variable!
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
50 pos=repname==classfamily
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
51 RC=unname(sum(lengths[pos])/totalLength)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
52 RC
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
53 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
54 repScore=function(repname){ # uses external variable!
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
55 pos=repname==classfamily
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
56 S=unname(mean(score[pos]))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
57 S
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
58 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
59
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
60 N=length(RMfiles)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
61
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
62 for (i in seq_along(RMfiles)){
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
63 cat(paste((RMfiles[i]),"\n"))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
64
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
65 sdir=getsubdir(RMfiles[i])
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
66 seqs.summary=fasta.summary(fasta[i])
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
67 totalLength=seqs.summary["length"]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
68
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
69 total.counts=data.frame(c("All_Reads_Length","All_Reads_Number"),c(NA,NA),c(NA,NA),c(totalLength,seqs.summary['N']),c(totalLength,seqs.summary['N']),c(NA,NA))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
70 colnames(total.counts)=c("class/fammily","class","family","hits","content","Mean_Score")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
71
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
72
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
73 # get RM data if not empty
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
74 rmsk_empty=length(grep("There were no repetitive sequences detected",scan(RMfiles[i],what=character(),n=1,sep="\n",quiet=T),fixed=T))>0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
75 rmsk_notperformed = length(grep("No Repbase used with RepetMasker",scan(RMfiles[i],what=character(),n=1,sep="\n",quiet=T),fixed=T))>0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
76 if (!rmsk_empty & !rmsk_notperformed) {
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
77 rmsk=read.table(RMfiles[i],header=F,as.is=T,skip=2,comment.char="*",fill=T)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
78 score=rmsk[,1]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
79 Ids=rmsk[,5]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
80 starts=rmsk[,6]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
81 ends=rmsk[,7]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
82 lengths=ends-starts
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
83 classfamily=rmsk[,11]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
84 out.table=table(classfamily)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
85 out.table=data.frame(names(out.table),c(out.table),stringsAsFactors=F)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
86 contents=sapply(out.table[,1],repCont)*100
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
87 meanScore=sapply(out.table[,1],repScore)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
88 class=gsub("/.*","",out.table[,1])
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
89 families=gsub(".*/","",out.table[,1])
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
90 out.table=data.frame(out.table[,1],class,families,out.table[,-1],contents,meanScore,stringsAsFactors=F)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
91 colnames(out.table)=c("class/fammily","class","family","hits","content","Mean_Score")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
92 out.table=out.table[order(out.table$content,decreasing=T),]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
93 out.table=rbind(out.table,total.counts)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
94 }else{
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
95 out.table=total.counts
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
96 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
97 out.table=out.table[order(out.table$content,decreasing=T),]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
98
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
99 write.table(out.table,file=paste(sdir,opt$output,sep=""),row.names=F,quote=F,sep="\t")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
100
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
101
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
102 # merge all tables:
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
103
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
104 colnames(out.table)=paste(colnames(out.table),c("",rep(sdir,5)))
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
105
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
106 if (i==1) {
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
107 mergedTable=out.table[,c(1,4,5,6)]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
108 }else{
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
109
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
110 mergedTable=merge(mergedTable,out.table[,c(1,4,5,6)],by=1,all=T)
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
111 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
112
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
113 }
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
114 mergedTable=mergedTable[order(rowSums(mergedTable[,seq(2,N*3,by=3),drop=FALSE],na.rm=T),decreasing=T),]
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
115 mergedTable[is.na(mergedTable)]<-0
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
116
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
117 write.table(t(mergedTable),file=paste(opt$output,"summary.csv",sep=""),col.names=F,row.names=T,quote=F,sep="\t")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
118
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
119 #save.image("tmp.RData")
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
120
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
121
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
122
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
123
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
124
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
125
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
126
a4cd8608ef6b Uploaded
petr-novak
parents:
diff changeset
127