Added functions winnow and winnow.mark to R package, updated manual

README.md

@@ -2,7 +2,7 @@
 
 Contact: Brandon Seah ([email protected])
 
-Cite: Brandon Seah (2014), genome-bin-tools, Online: https://github.com/kbseah/genome-bin-tools
+Cite: Brandon Seah (2014,2015), genome-bin-tools, Online: https://github.com/kbseah/genome-bin-tools
 
 ## 0. Introduction
 
@@ -36,6 +36,8 @@ Use your favorite assemblier, like IDBA-UD (http://i.cs.hku.hk/~alse/hkubrg/proj
 
 BBmap is convenient because it calculates length and GC along with the coverage per scaffold, and outputs them all to the same file. If you use a different mapping tool, you will have to aggregate all these data together in a tab-separated table with the same header names as the bbmap output.
 
+NOTE: Newer versions of BBmap (after 33.x) insert a comment character `#` in the header line of the covstats file. Remove this character before running the scripts.
+
 For differential coverage binning, you will need a second read library from a different sample where at least some of the genomes present in the first sample should also be present. Map the second read library onto the same assembly and save the coverage values as a different file:
 
 ```
@@ -144,7 +146,23 @@ If you see a cluster of scaffolds which you would like to save as a bin, you can
 
 `bin` is now an object of class genomestatsbin. Type the name of the bin object to see a summary of the bin. If you imported the marker, SSU, and/or tRNA data, a summary of how many of each are contained in the bin will be reported (this is useful if the marker genes are typical single-copy genes, for example).
 
-### 3c. Fishing for connected contigs using Fastg files (experimental)
+### 3c. Taking subsets of genomic bins
+
+If you want to get a subset of scaffolds based on GC, coverage, or length, use the function winnow()
+
+```R
+ > bin2 <- winnow(d,cov=c(200,Inf)) # Return the subset of contigs in d which have coverage above 200
+ > bin2 <- winnow(d,gc=c(0.2,0.5),cov=c(100,200),len=c(1000,Inf)) # Return the subset of scaffolds in d that have GC between 20-50%, coverage 100 to 200, and length above 1000
+```
+
+If you want to get a subset of scaffolds containing marker genes that belong to a particular taxon, use the function winnow.mark()
+
+```R
+ > bin3 <- winnow.mark(d,param="Class",value="Gammaproteobacteria") # Returns all scaffolds containing a marker gene whose value for "Class" is "Gammaproteobacteria
+```
+The functions winnow() and winnow.mark() work for both genomestats and diffcovstats objects
+
+### 3d. Fishing for connected contigs using Fastg files (experimental)
 
 Fastg files are generated by newer versions of the SPAdes assembler, and contain contig connectivity information generated during the assembly process. These can be useful, e.g. to "fish" scaffolds that are from the same genome but which were inadvertently left out of the interactively chosen bin.
 
@@ -154,8 +172,8 @@ For now you will have to manually edit the `genome_bin_tools.r` file to specify
 
 Fish out a new bin from an old bin using a Fastg file:
 ```R
- > bin2 <- fastg_fishing.genomestatsbin(d,bin1,"path/to/fastg/file.fastg")
- > bin2 <- fastg_fishing.genomestatsbin(d,bin1,"path/to/fastg/file.fastg", save=TRUE,file="bin2.scaffolds.list") # Save list of scaffolds in new bin to external file
+ > bin4 <- fastg_fishing.genomestatsbin(d,bin1,"path/to/fastg/file.fastg")
+ > bin4 <- fastg_fishing.genomestatsbin(d,bin1,"path/to/fastg/file.fastg", save=TRUE,file="bin2.scaffolds.list") # Save list of scaffolds in new bin to external file
 ```
 
 You can compare the two bins by plotting them overlaid:
@@ -207,3 +225,20 @@ Identical syntax to `fastg_fishing.genomestats`:
  > Bin2 <- fastg_fishing.genomestats(D,Bin1,"/path/to/file.fastg")
 ```
 
+## 5. Mapping and reassembly
+
+Once you have your final bin, either the shortlist of scaffolds in that bin was exported with the save parameter when that bin was defined, or you can use the native R function write():
+
+```R
+ > write(as.character(bin1$scaff$ID),file="shortlist.file") # for genomestatsbin object
+ > write(as.character(bin2$diffcov$ID),file="shortlist2.file") # for diffcovstatsbin object
+```
+
+This gives you a list of scaffold names in your bin. Use a tool like faSomeRecords to retrieve the corresponding Fasta sequence records. Use a mapper like BBmap to get reads that map to those scaffolds, and then reassemble with your favorite assembler.
+
+```
+ $ faSomeRecords original_assembly.fasta shortlist.file shortlist.fasta 
+ $ bbmap.sh ref=shortlist.fasta outputunmapped=f outm=reads_for_reassembly.fastq.gz in=original_reads.fastq
+```
+
+Good luck!

genome_bin_tools.r

@@ -2,7 +2,8 @@
 
 ## Tools for interactive genome binning in R
 
-## Version 1 - 2014-10-28 - Converted previous version to object-oriented paradigm
+## Version 1   - 2014-10-28 - Converted previous version to object-oriented paradigm
+## Version 1.2 - 2015-03-03 - Added functions winnow and winnow.mark
 ## Contact: [email protected]
 
 ## Required packages:
@@ -185,6 +186,37 @@ plot.genomestatsbin <- function(x, ... ) {              # inherit the same plot
     plot.genomestats (x, ...)
 }
 
+winnow <- function(x,gc=c(0,1),cov=c(0,Inf),cov2=c(0,Inf),len=c(0,Inf),save=FALSE,file="bin_scaffolds.list") {
+    # "Winnow" a genomestats(bin) or diffcovstats(bin) object by GC%
+    # Only return those contigs which are in that GC range
+    if (class(x) == "genomestatsbin" | class(x) =="genomestats") {
+        scafflist <- as.character(x$scaff$ID[which(x$scaff$Ref_GC >gc[1] & x$scaff$Ref_GC < gc[2] & x$scaff$Avg_fold>cov[1] & x$scaff$Avg_fold<cov[2] & x$scaff$Length>len[1] & x$scaff$Length<len[2])])
+        winnowedbin <- genomestatsbin (shortlist=scafflist, x=x, points=NA, save=save, file=file)
+        return(winnowedbin)
+    }
+    if (class(x) == "diffcovstatsbin" | class(x) == "diffcovstats") {
+        scafflist <- as.character(x$diffcov$ID[which(x$diffcov$Ref_GC > gc[1] & x$diffcov$Ref_GC < gc[2] & x$diffcov$Avg_fold_1>cov[1] & x$diffcov$Avg_fold_1<cov[2] & x$diffcov$Avg_fold_2>cov2[1] & x$diffcov$Avg_fold_2<cov2[2] & x$diffcov$Length>len[1] & x$diffcov$Length<len[2])])
+        winnowedbin <- diffcovstatsbin (shortlist=scafflist, x=x, points=NA, save=save, file=file)
+        return(winnowedbin)
+    }
+    else {cat ("Object must be of class genomestats, genomestatsbin, diffcovstats, or diffcovstatsbin!\n")}
+}
+
+winnow.mark <- function(x,param="Class",value="Gammaproteobacteria",save=FALSE,file="bin_scaffolds.list") {
+    # "Winnow a genomestats or diffcovstats by its marker table
+    # Only return those scaffolds whose marker parameter "param" is of value "value" (e.g. only return all scaffolds with markers whose Class designation is Gammaproteobacteria)
+    if (class(x) == "genomestatsbin" | class(x) == "genomestats") {
+        scafflist <- as.character(x$mark$scaffold[which (x$mark[,which(names(x$mark)==param)] == value)])
+        winnowedbin <- genomestatsbin (shortlist=scafflist, x=x, points=NA, save=save, file=file)
+        return(winnowedbin)
+    }
+    if (class(x) == "diffcovstatsbin" | class(x) == "diffcovstats") {
+        scafflist <- as.character(x$mark$scaffold[which (x$mark[,which(names(x$mark)==param)] == value)])
+        winnowedbin <- diffcovstatsbin (shortlist=scafflist, x=x, points=NA, save=save, file=file)
+        return(winnowedbin)
+    }
+    else {cat("Object must be of class genomestats, genomestatsbin, diffcovstats, or diffcovstatsbin!\n")}
+}
 
 reslicebin <- function(x=NA,bin,save=FALSE,file="bin_scaffolds.list") {
 # Reslice a bin so that it can be plotted on a diffcovstats or genomestats plot generated from an object different to the one originally used to produce the bin (get it?)