# set filenames for inout and output files blast.output.filename <- "~/Desktop/DD/m110923_025602_00123_c100181102555600000315046911131143_s2_p0.bas.h5_mapped_output.txt" output.filename <- "~/Desktop/DD/alignment_summary.txt" # read in blast alingment in tab delimited format the.data <- read.table(blast.output.filename, as.is=T, header=T) # generate summary count table for sseqid column count.table <- summary(as.factor(the.data[,2])) count.table <- cbind(BarcodeName=names(count.table), Count=as.numeric(count.table)) # initiate variables for basic information about sseqid matches med.length <- rep(NA, nrow(count.table)) max.length <- rep(NA, nrow(count.table)) percent.20 <- rep(NA, nrow(count.table)) barcode <- rep(NA, nrow(count.table)) # for each sseqid, generate summary information for (index in 1:nrow(count.table)) { temp.data <- which(the.data[,2]==count.table[index,1]) med.length[index] <- median(as.numeric(the.data[temp.data,"length"], na.rm=T)) max.length[index] <- max(as.numeric(the.data[temp.data,"length"]), na.rm=T) percent.20[index] <- length(which(as.numeric(the.data[temp.data,"length"])>=20))/length(temp.data) barcode[index] <- the.data[temp.data[(which(as.numeric(the.data[temp.data,"length"])==max.length[index]))[1]],"qseq"] } # generate final table and write output count.table <- cbind(count.table, med.length, max.length, percent.20, barcode) write.table(count.table, output.filename, sep="\t", row.names=F, col.names=T)