Database construction Filter1 step

[hemprichbennett]

Hi Devon, firstly thanks for a really interesting paper and making the tidybug code so open and well documented!

I'm trying to run the steps outlined in database_construction.md and am getting a bit confused. I've successfully used bold_datapull_arths.R to query BOLD. However then the Filter1 step starts by using grep on a file 'boldCustom.allArth.seqNtaxa.csv' to remove any lines matching 'sequenceID,taxon,nucleotides'. However none of the bold_datapull scripts in the repo output a file with that name, and I don't think any of them output a csv with column names matching that pattern? The only two csvs output with similar names are boldCustom.allArth.seqs.csv and boldCustom.allArth.taxa_string.csv, but using the command from Filter1 on either of these generates an empty fasta.

Thanks very much
Dave

[devonorourke]

Hammerhead twitter - glad to know someone else is out there in the world who cares about bat poop 😄

If you run bold_datapull_arths.R, it should create a .csv file called boldCustom.allArth.seqs.csv at this line. Likewise, you should get a file named boldCustom.allArth.taxa_string.csv a little later in the script here.

Clearly you don't get any file named boldCustom.allArth.seqNtaxa.csv, so I must have completely forgotten to include a step in the documentation... ooof. sorry 🤦

The boldCustom.allArth.seqNtaxa.csv is intended to be a concatenation/mergint of the information in those two .csv outputs: it'll include the sequenceID and nucleotides information boldCustom.allArth.seqs.csv, and merge it with the sequenceID and taxon information from boldCustom.allArth.taxa_string.csv.

As I wrote these scripts over 3 years ago, I honestly can't remember whether I merged these in R after the fact, or if I merged them with Bash. I think the simplest way to resolve this is to perform a merge within R and produce the necessary .csv file that fits with the documentation steps for filtering that data. I've just updated the bold_datapull_arths.R script to do that (see changes here) - mind testing that out and seeing if it works?

As you might imagine, I am not an experience programmer, and designing this database required me cobbling together a bunch of different bits of software. It's a major drawback, as I'd love to have enough Python chops to just make it all easily programmable for a user to enter whatever parameters they want and click "go", and have a database. Fortunately, the smarter folks behind QIIME are helping me do just that, so keep an eye out for an upcoming tutorial that is going to greatly simplify this COI database process in the near future - see here for starters on how it'll generally work (though it's for 16S in this example), and here for an ongoing discussion of how I'm doing this for COI.

[hemprichbennett]

#Ha, I thought I recognised your name from twitter! I can't remember what your handle on there is now for some reason though...

Ah great, I would be very interested in the upcoming tutorial you're making with the QIIME folks! Could I ask you to send me it when its out please, just in case I don't see it? This all sounds extremely useful for some work that I'll be doing in future.

Cheers for the updates to the script, cobbling together lots of different bits of bash, python and R sounds very familiar! I've never used data.table before so I'm not familiar with it's error messages. It required x.taxon and x.seqs to be transformed into data.table objects before running setkey, using
setkey(as.data.table(x.seqs), sequenceID) setkey(as.data.table(x.taxon), sequenceID).
However then when running line 187 it gave the error Error in [.data.frame(x.seqs, x.taxon, nomatch = 0) : unused argument (nomatch = 0)
However I think I've managed to do the same thing with seqNtaxa <- dplyr::inner_join(x.taxon,x.seqs). The Filter1 step now appears to work on the csv that the R script outputs, hopefully I can try out the rest of the pipeline soon :)

[devonorourke]

Thanks for testing Dave.
Yep, the error you were seeing with data.table was because I forgot to declare the two objects being merged as data.table objects. The inner join with dplyr, or just base R merge, should both do the trick. I've updated the script once more in hopes that this resolves the bug.
It's strange this was in the code in the first place, as the other two R scripts used to pull either all chordate records, or all the non-arthropod/non-chordate records both used the simple merge statement I ended up amending this script with! Hopefully you shouldn't have any trouble with those should you choose to use either of them.
Let me know if any future problems persist.

[hemprichbennett]

Hey Devon, I see from your conversation on the QIIME support conversation that you're making progress with the new workflow! I just thought I'd update you on how things were going with this one.

It looks like a change in the csvtk software means that the commands need to be altered from -f 4 to -f 2 to work, then it appears to work fine. Unless my doing that has changed something to a logical error, I didn't dig into it too thoroughly as I've been unable to install QIIME on my computers (OSX, Linux and Windows, sigh), so have had to give up once I get to the QIIME steps. As support for QIIME1 has ceased it looks like its basically uninstallable, which is annoying as it kills some of my old code!

As it seems like the QIIME2 workflow is coming along nicely I think I may give up at this point, good luck with it :)