Skip to content

Commit 0718b14

Browse files
add utils from TE_Aalb repo
1 parent 52ef843 commit 0718b14

File tree

1 file changed

+35
-0
lines changed

1 file changed

+35
-0
lines changed

utils/query_coverage.R

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
#!/usr/bin/Rscript
2+
3+
# https://github.com/clemgoub/TE-Aid/blob/master/consensus2genome.R
4+
# https://github.com/clemgoub/TE-Aid/blob/master/Run-c2g.R
5+
6+
Args = commandArgs()
7+
blast_file = Args[6]
8+
cons_len = as.numeric(Args[7])
9+
out_file = as.character(Args[8])
10+
11+
cons_coverage=function(blast_file=NULL, cons_len=NULL, out_file=NULL){
12+
blast = read.table(blast_file, sep="\t")
13+
new_rownames <- c(paste(blast$V2,"-",as.character(blast$V8),":",as.character(blast$V9), sep = ""))
14+
15+
#make the coverage matrix
16+
coverage = matrix(rep(0, length(blast$V1)*as.numeric(cons_len)), byrow = T, ncol = as.numeric(cons_len))
17+
for(i in 1:length(blast$V1)){
18+
if (blast$V7[i] <= cons_len) {
19+
coverage[i,]<-c(rep(0, blast$V6[i]-1),rep(1, abs(blast$V7[i]-blast$V6[i])+1), rep(0, as.numeric(cons_len)-blast$V7[i]))
20+
}
21+
}
22+
23+
rownames(coverage) <- new_rownames
24+
coverage<-colSums(coverage)
25+
write.table(coverage, file=out_file, col.names=F, row.names=F)
26+
27+
# number of full copies
28+
full=blast[abs(blast$V6-blast$V7) >= 0.9*as.numeric(cons_len),]
29+
cat(nrow(full))
30+
}
31+
32+
cons_coverage(blast_file = blast_file,
33+
cons_len = cons_len,
34+
out_file = out_file
35+
)

0 commit comments

Comments
 (0)