Paired-end R1/R2 reads not merging in ddRAD assembly #step1

[manhbio]

Hello all,

I got an issue in step 1. when I ran 7 ddRAD samples with paired-end reads (14 .fq.gz files in total) using iPyRAD v0.9.107.

During Step 1, it seems that the R1 and R2 reads are not being merged:

ipyrad [v0.9.107]
Interactive assembly and analysis of RAD-seq data
Step 1: Loading sorted fastq data to Samples
[####################] 100% 0:01:49 | loading reads
14 fastq files loaded to 14 Samples

Here is the format of my sample filenames (_R1 and _R2 present for each sample):

CJ200827-283_R1.fq.gz
CJ200827-283_R2.fq.gz
CJ200920-860_R1.fq.gz
CJ200920-860_R2.fq.gz
CJ200921-949_R1.fq.gz
CJ200921-949_R2.fq.gz
CJ201008-352_R1.fq.gz
CJ201008-352_R2.fq.gz
CJ201011-273_R1.fq.gz
CJ201011-273_R2.fq.gz
CJ230609-364_R1.fq.gz
CJ230609-364_R2.fq.gz
CJ240425-188_R1.fq.gz
CJ240425-188_R2.fq.gz

Could you advise why the merging might not be happening, or if there is any additional configuration required in params file for ddRAD paired-end reads?

Many thanks
Manh

[isaacovercast]

Can you share your params file? The datatype should be 'pairddrad', if you only put 'ddrad' then it will treat it as single-end.

[manhbio]

Thank you very much. I have reset the params file as your suggestion. then it's running well right now.

Could I ask another question? When working with a large dataset (around 300 samples) using the de novo method, step 3 (clustering/mapping) takes a very long time. I’m running it on a Linux server with 36 cores, but step 3 has already taken more than a month. So can we try selecting a subset of representative samples from all populations first, then using the resulting output file (loci file => change to Fasta file) as a reference for the full dataset? If so, would this approach affect the results, such as the SNP calling or downstream analyses? Thank you.

[isaacovercast]

@manhbio That is not a bad idea, and in fact it is the approach we are taking for ipyrad2, which is in active development. The fact that it is taking so long might be an indication that there is something else going on that might be worth looking into first. For example, 2 things that can cause slowdowns are high error rates in the raw data and/or low amounts of available RAM. The first issue is easier to problem solve. How long are the reads? Have you taken a look at a few of the samples with fastqc? This will tell you how much low quality data there is and it will give you information you can use to parameterize the trim_reads parameter in step 2. Take a look at fastqc and let me know how the data looks.

The reference method you suggest is a good idea, but lets make sure the data is looking good before we try something more complicated.

[manhbio]

Hi @isaacovercast, thank you.

Based on your information, it seems that there might be a small problem with my job.

1. RAM issue: I have 256GB DDR4, but I didn't explicitly set memory allocation in the script when submitting the job on Linux before.

2. Raw dataset quality: I have checked the raw reads using FastQC before, and overall the quality looks good. The reads are paired-end 150 bp, with high base quality across most positions. Only the last few bases (144-150 bp) show slightly lower quality, but it seems minor. Should I trim around 5 bp at the beginning and end? I will set params file then.

I would appreciate any advice you might have in this case, as step 3 seems to be taking too much of my analysis time. And I have attached the fastqc file below. Please take a look at this file if necessary.

Thank you for your help.

multiqc_report.html

[isaacovercast]

It wouldn't hurt to trim the last 10 bases off the reads, but I also don't think it would make a huge difference. It looks like there is some adapter contamination, but that should be handled automatically by step 2. How many raw reads are there per sample? I am not sure about the behavior of your system if you don't explicitly set the RAM, but if you are underprovisioned on RAM it will definitely go very slow.

[manhbio]

I have around 7-10M reads per sample

[isaacovercast]

That is a lot. Well, I will say as for your first idea, something that I have done is taken fq files for 3-5 representative samples from the *_edits directory (preserves filtering/trimming/adapter removal) and then used something like spades to build a pseudo-reference, then plug this pseudo-reference in the reference_sequence parameter. I have done this and it seems to work quite well, it's also much faster if you feel like trying that way.

[manhbio]

Yes, I am doing it that way right now. Thank you very much. I think the reason is that I have a large number of reads per sample and I am working with a big dataset (300 samples), so clustering in step 3 takes a long time. I also tried analyzing 7 samples before, and it took around 2 days to complete all 7 steps.

I will set the RAM in the script and hope it will help support this step. Thank you for your help.

[annavmello]

Hi all! I am having the same problem, but with half of data (177 samples), the step 3 never ends and i don't know what might be going on. I am also using ddRAD pair-end reads, but my reads are much smallers with 2-3 millions per sample. I'm working in a cluster and using 32 cores for this step. Phred score is also good > 35 and the step2 went really well removing adapters.

I also ran some tests with a smaller number of samples (28 samples) using R1 and R2, and it finished within 9 hours.

This is the params that i'm using:

ipyrad params file (v.0.9.107)
ipfusco ## [0] [assembly_name]: Assembly name. Used to name output directories for assembly steps
/cfs/klemming/projects/snic/naiss2024-6-261/AVAM/ipfusco ## [1] [project_dir]: Project dir (made in curdir if not present)
## [2] [raw_fastq_path]: Location of raw non-demultiplexed fastq files
## [3] [barcodes_path]: Location of barcodes file
/cfs/klemming/projects/snic/naiss2024-6-261//AVAM/links/*.fastq.gz ## [4] [sorted_fastq_path]: Location of demultiplexed/sorted fastq files
denovo ## [5] [assembly_method]: Assembly method (denovo, reference)
## [6] [reference_sequence]: Location of reference sequence file
pairddrad ## [7] [datatype]: Datatype (see docs): rad, gbs, ddrad, etc.
TGCAG, CGG ## [8] [restriction_overhang]: Restriction overhang (cut1,) or (cut1, cut2)
5 ## [9] [max_low_qual_bases]: Max low quality base calls (Q<20) in a read
33 ## [10] [phred_Qscore_offset]: phred Q score offset (33 is default and very standard)
6 ## [11] [mindepth_statistical]: Min depth for statistical base calling
6 ## [12] [mindepth_majrule]: Min depth for majority-rule base calling
10000 ## [13] [maxdepth]: Max cluster depth within samples
0.85 ## [14] [clust_threshold]: Clustering threshold for de novo assembly
0 ## [15] [max_barcode_mismatch]: Max number of allowable mismatches in barcodes
2 ## [16] [filter_adapters]: Filter for adapters/primers (1 or 2=stricter)
35 ## [17] [filter_min_trim_len]: Min length of reads after adapter trim
2 ## [18] [max_alleles_consens]: Max alleles per site in consensus sequences
0.05 ## [19] [max_Ns_consens]: Max N's (uncalled bases) in consensus
0.05 ## [20] [max_Hs_consens]: Max Hs (heterozygotes) in consensus
4 ## [21] [min_samples_locus]: Min # samples per locus for output
0.2 ## [22] [max_SNPs_locus]: Max # SNPs per locus
8 ## [23] [max_Indels_locus]: Max # of indels per locus
0.5 ## [24] [max_shared_Hs_locus]: Max # heterozygous sites per locus
0, 0, 0, 0 ## [25] [trim_reads]: Trim raw read edges (R1>, <R1, R2>, <R2) (see docs)
0, 0, 0, 0 ## [26] [trim_loci]: Trim locus edges (see docs) (R1>, <R1, R2>, <R2)
p, s, l, n, k, a, G, g, u, v, t, m ## [27] [output_formats]: Output formats (see docs)
## [28] [pop_assign_file]: Path to population assignment file
## [29] [reference_as_filter]: Reads mapped to this reference are removed in step 3

and this setup for the cluster

#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --mem-per-cpu=885M

also putting the cores in the ipyrad command line: ipyrad -p params_fusco.txt -s 3 -c 32

Does anyone have any advice? I'm really struggling with this, and any help would be really appreciated.

[isaacovercast]

Did you check your reads for quality with fastqc? If there are lots of low quality bases toward the tail end of reads then this can slow down clustering (you can trim the bad parts using trim_reads in step 2).

The other thing to consider is you are allocating 885Mb per core but we recommend 4Gb per core. Under-provisioning RAM in step 3 can also cause it to go very slow. You are normally better off with fewer cores and more ram per core if the ram on your host is limited. Let me know if that helps.

[annavmello]

Hi @isaacovercast,

Thank you for the quick response.

I did check the quality of the reads; here is the MultiQC report in case it is useful:

multiqc_report_R1R2.html

Regarding memory, 888 MB per core is the limit I can currently request on the cluster. I have contacted the cluster support to see if there is a way to increase this, but I don’t know yet if it will be possible.

If increasing per-core memory isn’t an option, I believe my only alternative is to wait, right?

[isaacovercast]

Is 888MB per core the max you can assign a core at all? Or is it 888MB is the max you can assign if you use 32 cores? For example, if you use only 16 cores can you assign 1.7GB per core? Or with 8 cores then 4GB? If so I would reduce the number of cores and increase the memory per core. If not then there's not much you can do and you'll just have to wait it out (though 888MB per core is an awfully low limit these days).