Skip to content
Snippets Groups Projects
Commit 3fa64749 authored by Shaman Narayanasamy's avatar Shaman Narayanasamy
Browse files

Fix bug in R plotting script

parent ec3b4852
No related branches found
No related tags found
No related merge requests found
......@@ -159,9 +159,23 @@ annot <- readGff3(annot_file, isRightOpen = FALSE)
# extract only interesting information
print("Processing gff3 annotation file")
annot.1 <- as.data.frame(cbind(as.character(annot@annotation$seq_name), annot@.Data,
str_split_fixed(str_split_fixed(getGffAttribute(annot, "inference"), ",", 2)[,2], ":", 3)[,2]))
annot.1 <- as.data.frame(
cbind(
as.character(annot@annotation$seq_name),
annot@.Data,
str_split_fixed(str_split_fixed(
getGffAttribute(annot, "inference"),
",", 2)[,2], ":", 3)[,2]
)
)
colnames(annot.1) <- c("contig", "start", "end", "database")
# Remove rows with no start and stop positions
annot.1 <- annot.1[-is.na(annot.1$start),]
annot.1 <- annot.1[-is.na(annot.1$end),]
# Replace non-annotated gens with NA
annot.1$database[annot.1$database==""]=NA
annot.1[,2:3] <- sapply(annot.1[2:3], as.character)
......@@ -172,6 +186,7 @@ print("Calculating no. of genes per contig")
annot.2 <- as.data.frame(table(annot.1$contig))
colnames(annot.2)[1:ncol(annot.2)] <- c("contig","no_of_genes")
save.image()
# calculate gene lengths
print("Calculating length of genes from gff3 file")
annot.1 <- as.data.frame(cbind(annot.1, annot.1$end - annot.1$start + 1))
......@@ -179,8 +194,16 @@ colnames(annot.1)[ncol(annot.1)] <- "gene_length"
# aggregate table and calculate total gene lengths within contig
print("Calculating coding density of contigs")
annot.2 <- cbind(annot.2, aggregate(gene_length~contig, data=annot.1, FUN=sum)[,2])
colnames(annot.2)[ncol(annot.2)] <- "total_gene_length"
# Create temporary table
total_gene_length <- aggregate(gene_length~contig, data=annot.1, FUN=sum)
colnames(total_gene_length)[ncol(total_gene_length)] <- "total_gene_length"
annot.2 <- merge(annot.2,
total_gene_length,
by="contig",
all.y=FALSE)
# remove temporary table
rm(total_gene_length)
# create annotation table
print("Creating annotation table")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment