1
$\begingroup$

I have raw paired-end RNA-seq reads for two novel eukaryotic species. Some background: the reads represent a copepod (arthropod) species each. The mRNA for each read set was obtained by extracting mRNA from multiple individuals. I am not exactly sure what sequencer was used (some Illumina short read machine as far as I can tell; I am not directly involved in this project).

I'd like to assemble these reads into transcriptomes de novo as no reference genomes are available as far as I know; my hope is to assemble with Trinity and then annotate with Dammit using its default databases (and check for completeness w/ BUSCO).

However, it looks like something went wrong during the sequencing runs (or during sample preparation) because the sequence GC content plots from FastQC look like this:

fastqc_gc_content_plot

Here the curve with the tall peak is one sample, and the more flatter curve is the other. There are two closely related curves per sample because FastQC apparently seems to have recognized each read (forward and reverse) as a separate sample.

I first presumed that this is just rRNA contamination (I believe the mRNA was prepared via poly-A selection), and I filtered the reads against all of sortmerna's default databases (which I believe include sequences from SILVA and Rfam among others). Then I passed the sortmerna'd read through FastP to filter out adapters, paired read correction, and low complexity filtering (alongside whatever else FastP does by default). I had ran FastQC before FastP to check if the rRNA filtering did anything and it hadn't, so I also asked FastP to filter by length (--length_required 151) since there were small amounts of reads with lengths shorter than this value which I thought could be the source of the additional peaks in the data.

However, even after sortmerna and FastP together, nothing has really changed:

sortmerna_fastp_fastqc_gc_content_plot

(There are only two traces/curves here instead of four because the paired reads were apparently recognized properly by FastQC now.)

My plan was to run Khmer to clean the reads a little bit more and then assemble them. I did do a quick assembly with Trinity (with the raw reads) and the BUSCO scores were pretty bad (only ~ 25% of the arthropod BUSCO genes present in the assemblies). I was hoping I'd be able to get Dammit to annotate whatever was left after the cleaning, despite the poor BUSCO prognosis.

Now I'm not entirely sure what to do if I can't clean the reads to begin with, since preventing nonsensical contigs would be absolutely important for a decent annotation.

So this leads me to my questions:

  • Are these data even salvageable?
  • Is there anything else I can do to clean the data?
  • Is this even contamination to begin with that I'm seeing here? (I am presuming some kind of contamination because both samples share that peak at ca. 18% GC content, among other shared features.) Would mapping all the reads against each other and discarding the shared sequences help? If yes, is there a tool I can use for this?
  • Any other thoughts in this regard?

I searched through the posts here and on Biostars as well as I could but I wasn't really able to find anything helpful in this regard (hence my own post).

I'm very new to all of this, so any inputs, tips, advice, and/or guidance would be most appreciated.

$\endgroup$
  • $\begingroup$ Please can you provide the other metrics given by fastQC? eg the quality per base? I didn't use the latest version by i am suprised by the y-axis in %. In addition, how did you run Trinity? on a merged set of reads or did you perform 1 assembly for 1 sample? $\endgroup$ – thomas duge de bernonville Jun 30 at 4:43
  • $\begingroup$ @thomasdugedebernonville here is a link to the fastqc files: github.com/vragh/fastqcfiles . As for Trinity, it was an assembly each for each sample (with that particular pair of reads). The samples are not from the same species. Also I actually went over the FastQC results themselves for the first time (I've been viewing them through MultiQC) and I noticed I have some overrepresented sequences that need to be removed (including some non-mainstream adapters, by the looks of it). $\endgroup$ – Dunois Jun 30 at 17:47
  • $\begingroup$ I am sorry but the links to fastqcfiles seems broken. Can you just run a fastp (without filtering according to length, or set a lower size threshold, eg 36 nt) followed by Trinity? You may also give a try to Bridger or BinPacker which are faster than Trinity without quality loss. $\endgroup$ – thomas duge de bernonville Jun 30 at 19:05
  • $\begingroup$ also, can you please tell me how you got the data. An academic platform? a private company? $\endgroup$ – thomas duge de bernonville Jun 30 at 19:10
  • $\begingroup$ @thomasdugedebernonville I apologize, I had accidentally made the repository private. It should be accessible now. I'll try assembling after FastP with defaults but I doubt it will help as FastP barely made a dent in the quality of the reads as mentioned in the post. The data is from an academic collaborator btw. $\endgroup$ – Dunois Jun 30 at 20:12
1
$\begingroup$

according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them. Try to run fastp with the following options:

fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

and run FastQc again to check to efficacity of the trimming.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hi, so I gave this a try, and there isn't much of an improvement (if any at all). I've uploaded the post-FastP FastQC results to the Github repo I linked to in my other comment. $\endgroup$ – Dunois Jul 1 at 10:57
  • $\begingroup$ sorry I can't see them on the repo... but can you provide the json or html output of fastp too? thanks! $\endgroup$ – thomas duge de bernonville Jul 1 at 11:18
  • $\begingroup$ They hadn't gotten uploaded somehow. All done now. There should be two zipped folders in the repo now--one contains the FastQC reports and the other the FastP reports. $\endgroup$ – Dunois Jul 1 at 12:27
  • $\begingroup$ ok thanks. Indeed, no strong improvement. I think you may provide custom fasta sequences to fastp using the --adapter_fasta argument. You may add sequences from Clontech Universal Primer and Clontech SMART CDS kits because contamination most likely come from that. Besides, do you have an idea about the number of genes in the genomes you are interested in? I saw you had less than 10M of reads and i wonder if it is sufficient for a complete de novo transcriptome assembly. $\endgroup$ – thomas duge de bernonville Jul 1 at 12:51
  • $\begingroup$ I'll try that. Would it make sense to remove all the overrepresented sequences? As for the number of genes, no clue honestly. One of the two species here (an unidentified Calanus) has a genome available for a related species (C. finmarchicus) from which we could probably get an estimate of the gene count but I don't think the other one has anything. I doubt the assembly will be complete in any sense of the word, but any transcriptome is better than no transcriptome, right? $\endgroup$ – Dunois Jul 1 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.