Too much information for Phyluce / SPAdes?

[jb23590]

Dear Phyluce Team,

I am hoping there may be an solution to issue we have come across. This may be an easy fix (I hope so) but we believe we may be feeding either Phyluce or SPAdes too much information for it to give accurate UCE counts. We have 95 samples sequenced at a proposed depth of 40M PE reads, and enriched with UCE baits from this paper
https://www.sciencedirect.com/science/article/pii/S1055790320302165?casa_token=9MG8aNwNNm4AAAAA:OUgk-MFQ3hhlznIn2NyzdbKOXmrylU_HYU5ThgxhzejTNFEcapjF0TBNPYLsEDjp_cVXGQCud0X6

Raw reads counts were on average 25M (min 8M, max 56M). After SPAdes, assembled contig numbers were 288 thousand (min 48T, max 1M). These are roughly as we expected. However, when running phyluce_assembly_match_contigs_to_probes and phyluce_assembly_get_match_counts, UCE counts are greatly lower than anticipated. UCE contigs are 780 on average (min 474, max 1190). The number of UCE loci removed for matching
multiple contigs averaged 389 (min 127, max 1259).

I also ran phyluce_assembly_match_contigs_to_probes with the min-identity and -coverage values set to 70, as here:

phyluce_assembly_match_contigs_to_probes
--contigs path-to-spades/contigs
--probes Cowman_etal_APPENDIX_C-hexa-v2-final-probes.fasta
--output Plate1ProbeMatches70
--min-identity 70
--min-coverage 70

I thought this would reduce the number of removed loci but the UCE counts are the same as when run with default values. We wonder if downsampling the 40M PE reads before running spades would also help, but have not tried this yet. At the moment, we are a little stuck as we had expected upwards of 1500 UCE contigs per sample.

Thanks in advance,
Jason

[brantfaircloth]

In short, these issues sound more related to the data you have collected than the routines in phyluce. That said:

You could attempt to downsample and see if that helps - it is possible, but downsampling having a significant (positive) effect on contig number also suggests something odd about your data like degraded DNA, some degree of contamination, etc. Contamination with other organisms would be more expected when working with corals than other taxa.

You could also try to manually assemble with spades, which will let you fiddle with various parameters, changing how spades assembles the data (e.g. using spades in metaspades mode versus single cell mode, changing coverage cutoffs to reduce assembly of contaminating contigs, etc).

The results you are seeing could also be due to the enrichment conditions used - if too strict and using baits that are divergent from targets, you might see lower total #s of UCE contigs enriched. That typically cannot be fixed on the software end - it would require adjustment on the lab end (lower annealing temp, change of buffers/protocols, etc.).

[jb23590]

Thank you Brant - we had indeed initially thought this was the result of error in the enrichment (we did not perform this ourselves) but hope this is not the case as we cannot rescue it if so.

Thanks for suggesting manual SPAdes. The current script is run in a parallel run and looks partially like this:

#######

phyluce_assembly_assemblo_spades
--conf $sampleid/shortsamples.$sampleid.conf
--output /tmp/singlerunlocal.$sampleid/
--memory 2500
--cores 40
>> singlerunlocal.$sampleid.out

#######

Would you suggest adding:

    --cov-cutoff 2
    --careful

?
I know these parameters are used frequently in other coral work but am unsure exactly why the value of 2 is given.
Thank you

[brantfaircloth]

I'm not entirely sure what might be best. If it was me, I'd just pull some data for a small # of individuals, assemble manually with spades (e.g. directly), then format as phyluce expects, and run searches for UCE contigs. What phyluce is running is the following spades command (adjusted to your # cores and memory):

spades.py --careful
--sc
--threads 12
--memory 62
--cov-cutoff 5
--pe1-1 read1
--pe1-2 read2
--pe1-s singleton
-o assembly

You could try auto for cov-cutoff:

spades.py --careful
--sc
--threads 12
--memory 62
--cov-cutoff auto
--pe1-1 read1
--pe1-2 read2
--pe1-s singleton
-o assembly

You could also just bump the cov-cutoff used for phyluce and just see what that does (or change to auto). You have to do this in the phyluce.config file. The default is already set to 5.

You might also try using metaspades, which would look like:

spades.py --careful
--meta
--threads 12
--memory 62
--cov-cutoff 5
--pe1-1 read1
--pe1-2 read2
--pe1-s singleton
-o assembly

Which turns on metaspades and sets cov-cutoff to 5.

You could even try mitofinder, which basically does a version of what's just above.

[jb23590]

I'll try these suggestions, thanks Brant. If you wondered why we are using parallel runs in the first place this is because when not in parallel (tried running both manual + reformatting, and through Phyluce) SPAdes was taking up so much cpu power that we were getting on average 0.5 samples a month on a node with 80 cores assigned.

To confirm, is the config file you refer to where we can change cov-cutoff the same where you list [samples] ?

[brantfaircloth]

You are welcome. And, that's totally bizarre re: CPU usage. Not sure of your setup, but things can sometimes get real weird for odd reasons. I usually run spades in parallel w/ 10-12 CPUs assigned per taxon and 4-6 of those jobs per node. We very often do it this way and don't even use our own wrapper within phyluce.

[jb23590]

Interesting - have you found a benefit of assigning the cutoff param manually rather than through your own wrapper? As the current script I have *using Phyluce) took weeks of adjustments to get it to work, i'd like to make minimal edits where possible.

Currently, it looks like this:

#######################################################################

Create main config file

cat -> Plate1Assembly.conf

[samples]
BH_BTHA_01:/home/jb23590/Plate1/Plate1Clean/BH_BTHA_01/split-adapter-quality-trimmed/
BH_BTHA_02:/home/jb23590/Plate1/Plate1Clean/BH_BTHA_02/split-adapter-quality-trimmed/
#etc

############################

mkdir Plate1SpadesParallel

awk -F : '{print "echo "[samples]" > ./" $1 "/shortsamples." $1 ".conf" }' Plate1Assembly.conf

Run output

awk -F : '{print "echo "" $1 ":/home/jb23590/Plate1/Plate1Spades/" $1 "/data/" >> ./" $1 "/shortsamples." $1 ".conf" }' Plate1Assembly.conf

Run output

awk -F : '{print "./singlerunlocal.sh " $1 }' Plate1Assembly.conf

Copy to text file for screen sessions later

############################

cat -> singlerunlocal.sh

#!/bin/bash

FILE: runSample.sh

sampleid="$1"
touch singlerunlocal.$sampleid.out
tnow=$(date)
echo "##### START TIME $tnow #####" >> singlerunlocal.$sampleid.out

source /usr/local/miniconda/bin/activate phyluce-1.7.1

phyluce_assembly_assemblo_spades
--conf $sampleid/shortsamples.$sampleid.conf
--output /tmp/singlerunlocal.$sampleid/
--memory 2500
--cores 40
>> singlerunlocal.$sampleid.out

mv /tmp/singlerunlocal.$sampleid/ ~/Plate1/Plate1Spades/Plate1SpadesParallel/$sampleid/
conda deactivate
tnow=$(date)
echo "##### END TIME $tnow #####" >> singlerunlocal.$sampleid.out
mail -s"**** SPADES ASSEMBLY JOB $sampleid COMPLETE"

#######################################################################

I will try first to add -careful and -cov-cutoff here before manually, but am not sure if this should be in the 'main config file' or --conf $sampleid/shortsamples.$sampleid.conf \

[brantfaircloth]

Not generally - cutoff of 5 works pretty good for us - get's rid of the crap and keeps most of the good stuff.