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$
3
  • $\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$ Jun 30, 2020 at 4:43
  • $\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$ Jun 30, 2020 at 19:05
  • $\begingroup$ also, can you please tell me how you got the data. An academic platform? a private company? $\endgroup$ Jun 30, 2020 at 19:10

1 Answer 1

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.

$\endgroup$
4
  • $\begingroup$ sorry I can't see them on the repo... but can you provide the json or html output of fastp too? thanks! $\endgroup$ Jul 1, 2020 at 11:18
  • $\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$ Jul 1, 2020 at 12:51
  • $\begingroup$ yes you're definitively! Concerning the overrepresented sequences, I never faced this problem for now. But the % is quite high (more than 3%) so i guess it deserved to be tested. As an alternative, you may try to map your RNA-seq reads on the related genome species (use HiSAT2+Stringtie pipeline). It might help to clean your reads! $\endgroup$ Jul 1, 2020 at 13:01
  • $\begingroup$ for sure, do not hesitate to contact me again. Hope you will find a solution for your data anyway. $\endgroup$ Jul 1, 2020 at 13:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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