massive 5.5TB uniq.k.2.c.2.seqs file - what's wrong?

[Melanie LaCava]

Hello,

I am trying to run dDocent-2.9.4 on a cluster computer with the RPE method, but after repeated attempts, I keep getting insanely huge output files for the second cutoff step in the process to reduce the dataset for cdhit. In my most recent attempt, you can see in the file list below that the uniq.k.2.c.2.seqs file got to 5.5TB(!!) before I canceled the job because clearly something is wrong.

The uniq.seqs file is quite large as well (93GB), and when I scrolled through it, I saw that although the counts column seemed accurate, seqs with more than 1 copy were still listed multiple times in uniq.seqs, when my understanding is that the seq should be listed only once, whatever the count. For example, one seq has a count of three, and in the uniq.seqs file, it is listed 3 times:

3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG
3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG
3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG

Shouldn't it only be listed one time? This applies to all multiples in that file, so it's not surprising the file is still large after the first data cutoff. Similarly, the uniq.k.2.c.2.seqs has lines repeated that I think should only be there once, for example, the following 4 lines are identical:

2 AACAAGAAAAAGCCAGCTGCGGCCGCAGCCGCCGCCGCCCCCGCGCCACCCTCTGACAGCGACAGCGACAAAGAGAGCGGCACCGAGAAAGACTCTGAGGACTCTGGCAGCTCGAACAAAGATCGGAAGAGCACACGTCTGNNNNNNNNNNTTGTTCGAGCTGCCAGAGTCCTCAGAGTCTTTCTCGGTGCCGCTCTCTTTGTCGCTGTCGCTGTCAGAGGGTGGCGCGGGGGCGGCGGCGGCTGCGGCCGCAGCTGGCTTTTTCTTGTTCTTGTCGTCCAGATCGGAAGAGCGTCGTGTAG
2 AACAAGAAAAAGCCAGCTGCGGCCGCAGCCGCCGCCGCCCCCGCGCCACCCTCTGACAGCGACAGCGACAAAGAGAGCGGCACCGAGAAAGACTCTGAGGACTCTGGCAGCTCGAACAAAGATCGGAAGAGCACACGTCTGNNNNNNNNNNTTGTTCGAGCTGCCAGAGTCCTCAGAGTCTTTCTCGGTGCCGCTCTCTTTGTCGCTGTCGCTGTCAGAGGGTGGCGCGGGGGCGGCGGCGGCTGCGGCCGCAGCTGGCTTTTTCTTGTTCTTGTCGTCCAGATCGGAAGAGCGTCGTGTAG
2 AACAAGAAAAAGCCAGCTGCGGCCGCAGCCGCCGCCGCCCCCGCGCCACCCTCTGACAGCGACAGCGACAAAGAGAGCGGCACCGAGAAAGACTCTGAGGACTCTGGCAGCTCGAACAAAGATCGGAAGAGCACACGTCTGNNNNNNNNNNTTGTTCGAGCTGCCAGAGTCCTCAGAGTCTTTCTCGGTGCCGCTCTCTTTGTCGCTGTCGCTGTCAGAGGGTGGCGCGGGGGCGGCGGCGGCTGCGGCCGCAGCTGGCTTTTTCTTGTTCTTGTCGTCCAGATCGGAAGAGCGTCGTGTAG
2 AACAAGAAAAAGCCAGCTGCGGCCGCAGCCGCCGCCGCCCCCGCGCCACCCTCTGACAGCGACAGCGACAAAGAGAGCGGCACCGAGAAAGACTCTGAGGACTCTGGCAGCTCGAACAAAGATCGGAAGAGCACACGTCTGNNNNNNNNNNTTGTTCGAGCTGCCAGAGTCCTCAGAGTCTTTCTCGGTGCCGCTCTCTTTGTCGCTGTCGCTGTCAGAGGGTGGCGCGGGGGCGGCGGCGGCTGCGGCCGCAGCTGGCTTTTTCTTGTTCTTGTCGTCCAGATCGGAAGAGCGTCGTGTAG

When I looked in the dDocent executable, I see that the RPE method has this special_uniq function that makes it hard for me to troubleshoot what might be the problem here. 

Any ideas for how to fix this? Requested info below.

Thanks!

Melanie

dDocent.runs

Variables used in dDocent (version 2.9.6) run at Fri Sep 22 11:04:22 PDT 2023
Number of Processors
8
Trimming
yes
Assembly?
yes
Type_of_Assembly
RPE
Clustering_Similarity%
0.9
Minimum within individaul coverage level to include a read for assembly (K1)
2
Minimum number of individuals a read must be present in to include for assembly (K2)
2
Mapping_Reads?
yes
Mapping_Match_Value
1
Mapping_MisMatch_Value
4
Mapping_GapOpen_Penalty
6
Calling_SNPs?
yes
Email

ls -l [NOTE: I was unable to paste the entire list (webpage became unresponsive) - there are many more samples with the same files and roughly the same file sizes as the two shown here]

-rw-rw-r-- 1 mlacava mlacava   57M Sep 22 11:19 70704_15.F.fq.gz
-rw-rw-r-- 1 mlacava mlacava   62M Sep 22 13:20 70704_15.R.fq.gz
-rw-rw-r-- 1 mlacava mlacava   61M Sep 22 16:52 70704_15.R1.fq.gz
-rw-rw-r-- 1 mlacava mlacava   67M Sep 22 16:52 70704_15.R2.fq.gz
-rw-rw-r-- 1 mlacava mlacava  1.4K Sep 22 16:52 70704_15.trim.log
-rw-rw-r-- 1 mlacava mlacava  218M Sep 22 12:09 70704_15.uniq.seqs
-rw-rw-r-- 1 mlacava mlacava   31M Sep 22 11:19 70704_16.F.fq.gz
-rw-rw-r-- 1 mlacava mlacava   34M Sep 22 13:20 70704_16.R.fq.gz
-rw-rw-r-- 1 mlacava mlacava   34M Sep 22 16:52 70704_16.R1.fq.gz
-rw-rw-r-- 1 mlacava mlacava   36M Sep 22 16:52 70704_16.R2.fq.gz
-rw-rw-r-- 1 mlacava mlacava  1.4K Sep 22 16:52 70704_16.trim.log
-rw-rw-r-- 1 mlacava mlacava  119M Sep 22 12:09 70704_16.uniq.seqs
-rw-rw-r-- 1 mlacava mlacava   417 Sep 15 16:14 config.file
-rw-rw-r-- 1 mlacava mlacava   478 Sep 22 11:04 dDocent.runs
-rw-rw-r-- 1 mlacava mlacava   54K Sep 22 15:11 lengths.txt
-rw-rw-r-- 1 mlacava mlacava   476 Sep 19 09:06 my_ddocent.sh
-rw-rw-r-- 1 mlacava mlacava  3.4K Sep 22 11:04 namelist
-rw-rw-r-- 1 mlacava mlacava  1.9M Sep 22 16:53 slurm-7525263.out
-rw-rw-r-- 1 mlacava mlacava  1.8M Sep 22 16:53 temp.LOG
-rw-rw-r-- 1 mlacava mlacava   27G Sep 22 16:07 total.f.uniq
-rw-rw-r-- 1 mlacava mlacava   72G Sep 22 15:18 total.fr
-rw-rw-r-- 1 mlacava mlacava   35G Sep 22 13:43 total.u.F
-rw-rw-r-- 1 mlacava mlacava   37G Sep 22 13:52 total.u.R
-rw-rw-r-- 1 mlacava mlacava   72G Sep 22 13:32 total.uniqs
-rw-rw-r-- 1 mlacava mlacava     0 Sep 22 11:04 trim.log
drwxrwsr-x 2 mlacava mlacava   770 Sep 22 16:52 trim_reports/
-rw-rw-r-- 1 mlacava mlacava  5.5T Sep 25 09:24 uniq.k.2.c.2.seqs
-rw-rw-r-- 1 mlacava mlacava   93G Sep 22 12:21 uniq.seqs
-rw-rw-r-- 1 mlacava mlacava   27G Sep 22 13:02 uniqCperindv

head temp.LOG [I have to cancel the job because the uniq.k.2.c.2.seqs file was getting out of control, so I just have the temp.LOG. This LOG has thousands of lines about perl locale error, so I cannot paste the whole log here. Here is the top up to the first perl error]

dDocent version 2.9.6 started Fri Sep 22 11:04:22 PDT 2023


At this point, all configuration information has been entered and dDocent may take several hours to run.
It is recommended that you move this script to a background operation and disable terminal input and output.
All data and logfiles will still be recorded.
To do this:
Press control and Z simultaneously
Type bg and press enter
Type disown -h and press enter
Now sit back, relax, and wait for your analysis to finish

Trimming reads and simultaneously assembling reference sequences
perl: warning: Setting locale failed.

[Melanie LaCava]

PS - I replicated the problem on a smaller subset of data. I took fastqs for just two individuals and ran through the data reduction steps manually (grabbing code from the dDocent executable file), and similarly ended up with a uniq.k.2.c.2.seqs file that is much larger than any of the intermediate files - shouldn't this file be the smallest since the sequences have been filtered twice? uniq.k.2.c.2.seqs. has 7469650 lines, but only 2640215 are unique - shouldn't this file only contain unique lines? This cascades to totaluniqseqs, and ultimately uniq.fq, which contains 7469650 sequences, even though those aren't all unique. Is that a problem?

ls -l

-rw-rw-r-- 1 mlacava mlacava  57M Sep 25 12:31 70704_15.F.fq.gz
-rw-rw-r-- 1 mlacava mlacava  62M Sep 25 12:31 70704_15.R.fq.gz
-rw-rw-r-- 1 mlacava mlacava  38M Sep 25 12:31 70704_15.f.uniq
-rw-rw-r-- 1 mlacava mlacava 107M Sep 25 12:31 70704_15.f.uniq.e
-rw-rw-r-- 1 mlacava mlacava 101M Sep 25 12:31 70704_15.forward
-rw-rw-r-- 1 mlacava mlacava 210M Sep 25 12:31 70704_15.fr
-rw-rw-r-- 1 mlacava mlacava 109M Sep 25 12:31 70704_15.r
-rw-rw-r-- 1 mlacava mlacava 109M Sep 25 12:31 70704_15.reverse
-rw-rw-r-- 1 mlacava mlacava 218M Sep 25 12:31 70704_15.uniq.seqs
-rw-rw-r-- 1 mlacava mlacava  31M Sep 25 12:06 70704_16.F.fq.gz
-rw-rw-r-- 1 mlacava mlacava  34M Sep 25 12:06 70704_16.R.fq.gz
-rw-rw-r-- 1 mlacava mlacava  22M Sep 25 12:31 70704_16.f.uniq
-rw-rw-r-- 1 mlacava mlacava  59M Sep 25 12:31 70704_16.f.uniq.e
-rw-rw-r-- 1 mlacava mlacava  55M Sep 25 12:31 70704_16.forward
-rw-rw-r-- 1 mlacava mlacava 114M Sep 25 12:31 70704_16.fr
-rw-rw-r-- 1 mlacava mlacava  59M Sep 25 12:31 70704_16.r
-rw-rw-r-- 1 mlacava mlacava  59M Sep 25 12:31 70704_16.reverse
-rw-rw-r-- 1 mlacava mlacava 119M Sep 25 12:31 70704_16.uniq.seqs
-rw-rw-r-- 1 mlacava mlacava   18 Sep 25 12:31 namelist
-rw-rw-r-- 1 mlacava mlacava 169M Sep 25 12:44 total.f.uniq
-rw-rw-r-- 1 mlacava mlacava 236M Sep 25 12:41 total.fr
-rw-rw-r-- 1 mlacava mlacava 114M Sep 25 12:41 total.u.F
-rw-rw-r-- 1 mlacava mlacava 122M Sep 25 12:41 total.u.R
-rw-rw-r-- 1 mlacava mlacava 236M Sep 25 12:40 total.uniqs
-rw-rw-r-- 1 mlacava mlacava 2.2G Sep 25 13:12 totaluniqseq
-rw-rw-r-- 1 mlacava mlacava 4.4G Sep 25 13:14 uniq.fq
-rw-rw-r-- 1 mlacava mlacava 2.3G Sep 25 13:13 uniq.full.fasta
-rw-rw-r-- 1 mlacava mlacava 2.2G Sep 25 12:48 uniq.k.2.c.2.seqs
-rw-rw-r-- 1 mlacava mlacava 336M Sep 25 12:31 uniq.seqs
-rw-rw-r-- 1 mlacava mlacava 169M Sep 25 12:32 uniqCperindv

[Jon Puritz]

Dear Melanie,

Believe it or not, those thousands of Perl warnings are important.  dDocent requires that LC_ALL=en_US.UTF-8.  I can't really help troubleshoot your cluster, but I would guess that if you are able to properly set the perl locale this will work properly for you.  

Here are some possible solutions: https://stackoverflow.com/questions/2499794/how-to-fix-a-locale-setting-warning-from-perl

Jon

--
You received this message because you are subscribed to the Google Groups "dDocent User Help Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ddocent+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ddocent/865104b3-4609-42cc-8f11-915d7031a56dn%40googlegroups.com.

--

Jon Puritz, PhD

Assistant Professor
Department of Biological Sciences
University of Rhode Island
120 Flagg Road, Kingston, RI 02881

Webpage: MarineEvoEco.com

Cell:    401-338-8739

Work:  401-874-9020

"The most valuable of all talents is that of never using two words when one will do.” -Thomas Jefferson

[Melanie LaCava]

Hi Jon,

Thanks for responding. I was able to get the perl locale settings fixed on my university cluster computer, and I no longer see errors for this in the temp.LOG. After this fix, I re-ran dDocent, and the same problem occurred - big 100GB uniq.seqs file with lines repeated, massive uniq.k.2.c.2.seqs file with many identical lines repeated, and I had to stop the job once it started to reach multiple TB in size.

When I went through the code from the dDocent executable manually, I noticed 2 unusual things:

Line 894:

AWK4='{for(i=0;i<$1;i++)print}'

This line seems to cause the uniq.seqs file to end up listing each sequence the number of times it is found, rather than listing it once with the count of occurrences in the first column. For example, I would expect the follow line to only be listed once, but instead it is listed 3 times, and I think that is due to this AWK4 expression - is this behaving correctly?

3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG
3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG
3 AGCACTCGATTTACATTTTCGCCGGCGGTGCACATCATTAAAAAGGCGGCGTGGCGACGCCTCATTAATCCCCTTCTGCGGCGACAATGGCGCCGCAAGAAGCGATCGACTCCATTCAAAAGC
GCGTGTCACGCTTTTGCCNNNNNNNNNNGTCGCGACCGTTCGCGGGTGCACAACTGTTCTTTTATTGCGGTCGCGACCGGAACCTCTCCTCTCTCGCCTCTCAACTTGGTGTTCCTCGCATGCATCGAAGC
AAAAGCGTGACACGCGCTTTTGAATGGAGTCGATCGCTTCTTGCGGCG

Line 993

join -1 2 -2 1 -o 1.1,1.2,2.2 total.f.uniq total.fr | mawk '{print $1 "\t" $2 "NNNNNNNNNN" $3}' | mawk -v x=$CUTOFF2 '$1 >= x' > uniq.k.$CUTOFF.c.$CUTOFF2.seqs

The input files here behave as I would expect in that total.f.uniq is identical to uniqCperindv (no repeated seqs), and these two files are a reasonable size (27GB and 72 GB, respecitvely). But the output of this line of code gets massive when creating the uniq.k.$CUTOFF.c.$CUTOFF2.seqs file.

Any thoughts?

Melanie

[Jon Puritz]

What type of RADseq are you analyzing?

To view this discussion on the web visit https://groups.google.com/d/msgid/ddocent/a6ca997e-ebc0-4a98-a003-aa7161556915n%40googlegroups.com.

--

Jon Puritz, PhD

Associate Professor

Department of Biological Sciences
University of Rhode Island
120 Flagg Road, Kingston, RI 02881

[Melanie LaCava]

The data were generated using BestRad with a Sbf1 enzyme and random shearing, so I selected RPE for my assembly method

[Jon Puritz]

Ok, before we go on a deep dive here.  Can you make sure to delete all *.seqs file from your working directory and restart?  It's possible that files from before (when the perl local wasn't right) have stuck around and are messing things up.  Let's make sure that's not the case first.

Jon

To view this discussion on the web visit https://groups.google.com/d/msgid/ddocent/fc895459-a29b-47e2-97a8-c027f41e7a4cn%40googlegroups.com.

[Melanie LaCava]

That looks to have fixed the problem, thanks!