Question

Visualization tool for GC-bias in SE data?

0

Entering edit mode

AKSR ▴ 20

@aksr-5026

Last seen 4.4 years ago

Need your help with 2 related questions, please. Thanks, in advance!

QUESTION 1

From this ALPINE link, under the section Description:

It is currently designed for un-stranded pairedend RNA-seq data.

So for both visualization and correction of GC-bias in RNA-Seq, I understand ALPINE is a preferred tool, but requires PE data. Right? Is there a GC-bias visualization tool for Single End Illumina RNA-Seq reads?

QUESTION 2

Is there a GC-bias correction tool for Single End Illumina RNA-Seq reads?

GC-bias ALPINE visualization • 1.8k views

ADD COMMENT • link updated 4.4 years ago by Michael Love 42k • written 4.4 years ago by AKSR ▴ 20

0

Entering edit mode

based on: https://www.biostars.org/p/446378/

ADD REPLY • link 4.4 years ago ATpoint ★ 4.5k

score 3 · Accepted Answer · 2020-06-28

3

Entering edit mode

Michael Love 42k

@mikelove

Last seen 1 day ago

United States

Alpine requires PE. For SE, you don’t know the size of the fragment, but can try to infer it. An approximation with longer SE reads is to just look at read GC. I always recommend MultiQC with the FASTQC GC curves and even wrote a plugin to draw a theoretical GC curve for human or mouse transcriptomes.

https://multiqc.info/docs/#theoretical-gc-content

ADD COMMENT • link 4.4 years ago Michael Love 42k

0

Entering edit mode

Thank you, this is helpful. A naive follow-up question, sorry:

What R syntax lines / package(s) do I use to directly input transcript FASTA file to generateDistn of your fastqcTheoreticalGC BioConductor package, to generate the final 2 column tsv file of %GC Vs. % of genome for use with MultiQC?

Does the syntax below look OK to you (I did generate a text file that looked OK, but still, to be 100% sure), THANKS!:

library("BSgenome")
library("fastqcTheoreticalGC")

Input_FASTA <- readDNAStringSet("Alltranscripts.fasta", "fasta")
generateDistn(seqs=Input_FASTA,
              file="fastqc_theoretical_gc.txt",
              name="MyFaveSpecies")

ADD REPLY • link 4.4 years ago AKSR ▴ 20

0

Entering edit mode

What organism are you using? If it's human or mouse, you don't have to generate the distribution, just use the built-in curves.

ADD REPLY • link 4.4 years ago Michael Love 42k

0

Entering edit mode

Medicago truncatula, a legume plant species

ADD REPLY • link 4.4 years ago AKSR ▴ 20

0

Entering edit mode

Ok so the first argument of generateDistn is actually described in the man page, see:

?generateDistn

It says you should input a DNAStringSet of transcript sequences. There are many places to obtain a transcript FASTA file and then you can read this into R using the Biostrings package to obtain something to provide to the first argument of generateDistn.

ADD REPLY • link 4.4 years ago Michael Love 42k

0

Entering edit mode

Yes, thank you. Since my original comment, I'd updated it to reflect the syntax I used ( I think successfully) to obtain the end result. I am not sure if the syntax lines using readDNAStringSet, followed by your generateDistn are visible to you. If they are visible, I request you please confirm if it looks OK to you, is what I am requesting you. Thanks, in advance.

ADD REPLY • link 4.4 years ago AKSR ▴ 20

1

Entering edit mode

Yes, this looks correct

ADD REPLY • link 4.4 years ago Michael Love 42k

0

Entering edit mode

Thanks a lot, I appreciate it.

ADD REPLY • link 4.4 years ago AKSR ▴ 20

0

Entering edit mode

With your help I added this dashed GC theoretical curve for the transcriptome for the species of interest. Thanks Michael - the green curves are for the ~ 150 libraries I want to check GC-bias for (single ended, cuz sorta old data)

In this image, does the shift of actual data (green) to the right of the theoretical GC curve (black dotted) mean that there is indeed some GC bias, resulting on-average in the preferential sequencing of genomic regions with a slightly higher GC % than would be expected if there were no bias? Could you please confirm or correct my interpretation? Thanks, in advance!

ADD REPLY • link 4.4 years ago AKSR ▴ 20

1

Entering edit mode

My guess is that there is probably not technical bias (which looks a little different, see alpine paper), but instead you are expressing some genes with higher GC content than what is computed by the theoretical curve (all transcripts at a constant rate). You can actually modify the "weights" per transcript in the script if you want to explore that - you could put in the estimated expression of each transcript instead of constant.

ADD REPLY • link 4.4 years ago Michael Love 42k