Hi everyone,
I am a PhD student and would greatly appreciate your help. I apologize in advance for my lack of knowledge.
I recently received a DNA sequencing data (llumina NovaSeq X sequençing, GBS or Genotyping-by-Sequencing libraries , enzyme PstI , paired-end 150 bp) for genetic diversity analysis and GWAS. As shown in the figure B below, the FastQC results revealed an abnormal GC-content distribution in the R2 reads, with a high proportion of poly-G reads (up to 33%). And this does not appear to be the typical poly-G tail issue; rather, many of the reads seem entirely of Gs.
I compared these data with a previous dataset generated in 2023 (Fig A) from the same species, using a similar GBS approach (Genotyping-by-Sequencing, PstI, paired-end 150 bp) but sequenced on Illumina NovaSeq 6000. This problem is not observed in that dataset.
I contacted the platform who did the library preparation and sequencing, and they suggested processing the data as single-end. Following this recommendation, I cleaned R1 and R2 separately using Cutadapt (to remove poly-G reads and adapters), concatenated them after trimming, mapped them to the reference genome using BWA, and performed SNP calling with STACKS. Unfortunately, the proportion of loci shared among individuals was very low (around 20%). After applying filters such as SNP missingness < 75% and MAF > 2%, I retained only about 2,500 SNPs, with an overall missing data rate still as high as 69% and samples missingness between 30% and 90% (although we had expected to find at least 40,000 SNPs after applying these filters).
I then tested several alternative pipelines (fastGBS, GATK, STACKS denovo pipeline) and even tried using only the R1 reads, since they do not contain these abnormal G-only reads, but I obtained the same results. In contrast, the older GBS dataset generated on the NovaSeq 6000 were doing very well.
Has anyone encountered a similar issue with NovaSeq X GBS data or RADseq data? Could this problem be related to the sequencing platform, library preparation, or some other technical artifact?
Any suggestions to process and save this data or feedback would be greatly appreciated.
Thank you very much for your help.
•
written
by
Butterfly •
0 answers
No answers yet.
Log in to answer this question.
-
Poly-G Poly-C issues and Adapter Trimming
written by Jon •
-
Why does including BQ<10 bases change Mutect2 calls in NovaSeq X Plus ctDNA data?
written by Anthony •
-
NovaSeq Xplus top unknown barcodes
written by 1769mkc
-
Tool: New Illumina error mode, new BBTools release (39.09) to deal with it
written by Brian Bushnell
-
Phylogeny from mixed GBS and WGS
written by maxrwjones •
-
Small RNA sequencing using Illumina 2 channel SBS: how to deal with Gs?
written by antoinefelden •
-
ATAC-seq troubleshoot - Just Noise
written by vk •
-
Kmer content failed in FastQC analysis
written by josmantorres •
-
what cause poly-G from NextSeq
written by CY
-
News: NovaSeq from Illumina (updated info)
written by GenoMax
Concatenated as in R1+R2 (joined the two files start-to-end) to treat as a single input? Are you certain one can use this as an input to STACKS (since the radtags would not be recognized properly).
33% of reads being pure poly-G sequence seems to be indicative of a problem in one (more or all) of the things above.
Have you tried removing these poly-G reads (and their corresponding counterparts from R1/R2 file) and then using the remainder of the data as paired-end input.
Thank you for your help
As far as I know , for GBS libraries, the PstI restriction sites do not have a fixed orientation. As a result, the same locus can be found in either R1 or R2. So by using only R1 reads , I will lose some informations from R2 reads. Therefore, discarding entire paired reads whenever R2 contained these poly-G reads may be alse overly conservative, especially if a large proportion of R1 reads are still usable. That's why I concatenated the cleaned R1 and R2 reads end-to-end in order to use them as a single-end input to retain as many informative reads as possible.
But Yes, I have already tried removing the poly-G reads together with their corresponding mates (with cutdapt) and then rerunning the analysis using the remaining reads as paired-end data. Unfortunately, this did not resolve the issue: the proportion of loci shared between samples remained very low, and the overall level of missing data was still very high.
Sounds like you have tried logical options. That coupled with above may indicate that this observation is a characteristic of your sample/libraries.
There are RADseq datasets in SRA done on NovaSeqX so it must work in general for reduced representation genome sequencing.
Since I obtained essentially the same results after trying several different approaches, I agree with you that the issue may be specific to our samples or libraries rather than to the downstream analysis alone. There are also recent papers reporting successful use of NovaSeq X for RADseq and GBS libraries.
So, I would mainly like to know whether anyone else has encountered a similar problem as me with NovaSeq X when working with RADseq or GBS libraries (a high proportion of poly-G reads and especially a high amount of missing data after variant calling), or I am the only one. It will be difficult for me to convince the reviewers with my data especially with this missing rate so I'm trying to find a solution until today.
I found some clues in Illumina’s documentation like Loading concentration considerations when migrating Illumina libraries between sequencing platforms, Library Loading Concentration and Optimization Guide for NovaSeq X/X Plus Instruments, and PhiX Spike-In Requirements for Low-Diversity Libraries on NovaSeq X Series Instruments. These documents suggest that loading concentration and low-diversity library (like GBS libraries) handling may be important factors . However, since I was not the person who performed the sequencing, I do not have access to the run parameters needed to evaluate these possibilities directly.
Do you need to compare the present data with one done on NovaSeq 6000? If that is the case then you have a bit of dilemma. You may end up having to redo the present libraries (if you are not 100% sure they are of good quality), optionally resequence them on NovaSeq 6000 or resequence old libraries on NovaSeq X (to evaluate platform "batch" effect, it should not exist, in theory). Unfortunately all of these options have costs (time and money) associated with them.
Had you told the sequence provider that these were GBS libraries? With reduced representation sequencing (low nucleotide diversity) loading conc can indeed be important. Sequence providers tend to pool multiple user samples together (since the throughput with platforms such as NovaX is so high and they need to achieve a price point for end-users). This can potentially lead to suboptimal results for samples that are not plain genomic sequence.
HI, I have been dealing with somewhat a similar but also different situation.
I also have reads done on the Novaseq 6000, LOTS of Poly-G issues. Some reads are all GGGG's, some just tails. A lot of the Poly-G is after the adapter, so those should get removed with adapter trimming. Then I also have entire reads that are Poly-G or Poly-C, I at first just tried removing these, but then realized there was another issue.
In my samples, I have three species that have long stretches of Poly-G or Poly-C within the reference sequences, so when I was trying to remove what I thought was artifact, was in fact my species. I also had interest if I was removing my sequences with my Adapter removal, so I chopped my reference fasta to 125bp and found out that I was removing a considerable amount of my sequences with the adapter removal.
I actually did one more thing that seemed to improve my mapping speed and accuracy. In all of my read sets, I ended up running into unpaired read issues that seemed strange. So I did an initial pre-processing step. I used FASTQtoSAM to convert my fastq into an unaligned BAM. Then using Samtools Fastx, extracting FASTQ from the BAM file, deinterlaced and back to a paired set again. (someone may decide this was wrong, but seemed to fix my issue)
So back to the first issue.
Since generally speaking, not all reads will have adapters attached, those are theoretically clean. I set both Hisat2 and Bowtie2 up for Perfect Mapping / 99.5% ID Mapping / 98% Mapping. First mapping my reference, then my host through Perfect/99%/98%, so this whittled down my original reads considerably. So at 98% there should at most be only a couple of bases of adapter present. This is where I did my adapter removal, followed by trimming ends to Q3 to remove any N's, and minimum length. Following this using BBmerge to error correct the R1/R2 and including the adapter to BBmerge, error correct without merging.
Then following with mapping again to 98%. This is where I removed the Poly-G using BBduk with K-29
I know it's a lot more work than should be necessary, but I think this finally solved this issue for me. Essentially separating the known good reads from the reads with issues so that any fix doesn't affect the good reads.
Moving this to a comment since the original question is about a specific technique called GBS.
Back to your post above:
One should never need to do extended manipulation of data to do relatively straight forward analysis. Since we can't see/access your data, we go on your assurance that what you did seemed to work in the end. When you get around to publishing the data, hopefully you will be able to convince the reviewers.
Hi GenoMax, here is an example.
The issue regarding Fastq to SAM. FASTQC showed all reads at 35-149bp for each read. I didn't notice a problem until I was mapping with BBmap, and it was leaving me with a lot of single reads. I'm assuming there is some format somewhere that lost the read pairing. I just found after converting to BAM and back to FASTQ that it fixed the problem, though not sure what the issue was. It technically shouldn't change any information, but only make sure each read is paired rather than losing it.
Taking the reference genomes that are known in my samples, I split them into 125bp. I did adapter removal with Trim Galore with standard defaults. The adapter removal is removing a significant amount of my dna, I had suspected this previously, but didn't realize how bad. This is with both Illumina Universal and Nextera. I'm currently just trying to map the reads to consecutively lower match %, then I should be able to use Tadpole to error correct any adapter in the tails.
I agree this shouldn't be that difficult, maybe you have a different workaround?
As I recall you are doing all this via
Galaxyso there is that additional layer. Let us not highjack this thread but what you describe above does not sound "normal".If you intend to publish the data eventually, keep in mind that your analysis procedure needs to be clear enough, so someone else could potentially replicate the result independently.
Yes, I'm not sure the issue. It's on Galaxy, but the tools are common. My assumption is that it does some sort of pairing Repair, similar to BBtools Repair.sh where there is some synchronization issues between some pairs.
Now with the image I added, these are ONLY reference genomes cut into 125bp pieces, I think you would agree that the Illumina and Nextera adapters have too much similarity to my references. I think this is where I ran into a different issue, if I had say 100,000 reads that map perfectly to reference, and I use BBmerge and include the adapter sequence as recommended, a good portion of reads aren't going to map perfectly anymore when BBmerge tries to correct the part of the dna that matches the adapter. So I think if I skip adding the known adapter into BBmerge that should fix that problem. Essentially ONLY using BBmerge for the adapter removal as standard adapter removal was butchering the reads.
I believe I found the issue with this, and it was due to using BBmap. I found a reference stating that BBmap is looking for /1, which is different than what is in the standard RAW reads.
@AV230601:CIDL1-0578:2339689248:1:10102:0394:0074 1:N:0:ATGGTTGTTG+CCTCATAACT
After converting FASTQ to BAM, back to FASTQ and de-interlaced
@AV230601:CIDL1-0578:2339689248:1:10102:0394:0074 /1
So BBmap was dropping what it believed was an unpaired read. So my workaround was to just separate my Host Reads from Other Reads, then convert the Other Reads to /1 /2, so that I don't have any issues with the other tools.
/1nomenclature was used a long time ago at the dawn of NGS for Illumina. If a read does not contain the "standard" Illumina ending for header that is common now, BBMap (or many programs that were around at the beginning) will look for/1or/2as a fall back, which is what your conversion from FASTQ --> BAM --> FASTQ did. Not sure how you lost the standard endings in Galaxy manipulations. Let us not discuss that here since it is unrelated to the original question.