1
$\begingroup$

Suppose I have single-end RNA-seq data for which the reads in the fastq file are reversed with respect to the original extracted RNAs.

Suppose I have the following workflow:

  1. Map the reads on the reference genome, sort and index the resulting sam output.
  2. Use the bam file to make a bigwig file showing reads mapping on the plus and minus strands of the chromosomes.
  3. Use the bam file and an annotation file to quantify the gene expression.

Steps 2 and 3 depend on step 1, but are not dependant on one another.

I'm sure that strandedness needs to be taken into account no later than at step 3.

The questions are the following:

  • Is it common to have mapping tools with an option to "correct" for strandedness?

  • Assuming the bam file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Regarding this second point, I'm split between two attitudes:

  1. There should be no correction, the bigwig files are there to represent the "raw" mapping results.
  2. The bigwig files are there to give an idea of the origin and abundance of the originally extracted RNAs, so a strandedness correction should be used if necessary.

Which one is the most common attitude?

If the second one is preferred, then wouldn't it be "safer" / "cleaner" to handle strandedness as early as possible? Possibly at the mapping, or even using a pre-processing step?

$\endgroup$
2
$\begingroup$

I usually take into account library strandness during the mapping to the reference step (but I never worked with bigwig files myself).

Assuming RNA-seq data from Illumina, you can use Hisat2 for alignment to the reference genome, and you have the option --rna-strandness to specify the strand-specificity information (default unstranded). This is particularly important when you wish to have XS strand attribute in your bam file for easier compatibility with downstream software (especially important in my experience if you are planning to use some utility from Cufflinks for your step 3).

$\endgroup$
4
  • $\begingroup$ Good to know that histat2 has such an option, but if I understand correctly, the aligned reads are not "corrected for strandedness": it is up to the next steps to do this, taking into account the XS attribute. Am I right? $\endgroup$
    – bli
    Jul 18 '17 at 18:02
  • 1
    $\begingroup$ @bli What would "corrected for strandedness" even mean? The aligners just use these options to preferentially produce results in a given manner (i.e., in a given orientation relative to genes). $\endgroup$
    – Devon Ryan
    Jul 18 '17 at 18:06
  • $\begingroup$ @DevonRyan I mean that, knowing that the reads are produced by a protocol giving the reverse of the RNAs, the aligner reverse them and report the mapping of the reversed reads. This would amount to pre-processing the reads, and then assume that the reads correspond to the RNAs in the rest of the pipeline. When you say "in a given orientation relative to genes", I suppose that implicitly, an annotation file is given to the mapper, and not just a genome index. $\endgroup$
    – bli
    Jul 18 '17 at 18:17
  • 1
    $\begingroup$ @bli Correct, the strandedness settings are only useful if the aligner is aware of the transcriptome (either due to directly aligning against it or having been given an annotation). $\endgroup$
    – Devon Ryan
    Jul 18 '17 at 18:21
2
$\begingroup$
  1. If your aligner has a strandedness option then go ahead and use it. The general idea here being that you'll preferentially align correctly to genes.
  2. Whether to bother with strandedness here depends on your goal. Are you interested in anti-sense transcription? Do you need to use the bigWig track to accurately quantify translational pausing or some other strand-specific phenomenon? If yes, then go ahead and make different bigWig files for each strand. If you just need the bigWig file for convenience in IGV or creating heatmaps (e.g., with deepTools) then you'll generally be fine without accounting for strand. For well-annotated species, the rate of anti-sense transcription is generally low and the rate of genes overlapping isn't absurdly high. Consequently, you can usually get away with using an unstranded bigWig.
  3. If you're using featureCounts you'll certainly want to use the appropriate strandedness setting.

It's relatively common for RNAseq aligners/mappers to offer strandedness options. Certainly tophat2 and hisat(2) do, but also things like salmon.

Again, for bigWig files the question becomes whether you really need two files to get the information you need. You're not going to delete the BAM files, so you can always make them if you need them.

$\endgroup$
1
$\begingroup$

Firstly, right at the start if the experiment it's important that your RNA samples are processed with a strand-specific protocol (e.g. Illumina's TruSeq Stranded) in order to produce stranded libraries for sequencing. If the samples haven't been treated as such, then there's nothing you can do to correct it.

Secondly, I'm not sure stranded protocols are compatible with single-end reads. The aligners I'm aware of use the direction of aligned read pairs to assess the strand from which the RNA fragment has come from. How that works with single-end reads, I don't know?

Finally, to answer your specific question, essentially step 1 is where you need to start with strand-specific options and do so at all subsequent steps as appropriate.

$\endgroup$
3
  • $\begingroup$ All stranded protocols are compatible with single-end reads. The library prep is identical, you just aren't flipping things and sequencing the other end of the fragment. $\endgroup$
    – Devon Ryan
    Jul 18 '17 at 17:46
  • $\begingroup$ @DevonRyan and at the feature counting stage? You can simply trust the reads are aligned to the correct strand, right? $\endgroup$
    – ithinkiam
    Aug 18 '17 at 8:48
  • 1
    $\begingroup$ You can trust that reads are aligned correctly given some MAPQ threshold, so if they have an orientation compatible with a given gene then yes you might as well believe that. $\endgroup$
    – Devon Ryan
    Aug 18 '17 at 8:50
1
$\begingroup$

Strandedness should be taken into account at the mapping stage.

If you are presented with only a BAM file, ask for the raw reads. If that's not possible, you can recover the reads (as fastq sequences) from the BAM file and/or remap them, but that conversion requires a lot more effort (e.g. most mappers expect the fastq files from each end to be sorted by read name and in the same order).

Regarding your specific questions...

Is it common to have mapping tools with an option to "correct" for strandedness?

Yes, mapping tools can be provided with strand direction information for mapping. The way that this is done depends on the program, for example bowtie2 has options --fr, --rf, and --ff to take strandedness into account. I haven't done any of my own comparisons to determine whether these options actually make a difference in mapping for Bowtie2. I know that Salmon (which does mapping and quantification at the same time) cares a lot more about strandedness, and will complain/warn if a stranded library type is discovered and not specified.

Assuming the BAM file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Yes, but it'd be better to remap with [the correct] strandedness information included. When I generate bigwig files to demonstrate strandedness (even assuming that strandedness is incorporated into the mapping), I still need to do a fair amount of legwork to convert that into bigwig files. The data needs to first be split into four groups (e.g. with samtools view -F 0xXX -f 0xXX):

  1. First read, mapped in the reverse orientation
  2. Second read, mapped in the reverse orientation
  3. First read, mapped in the forward orientation
  4. Second read, mapped in the forward orientation

I then combine two of those groups, depending on the library strandedness scheme. For a typical stranded Illumina run, groups 1 and 4 are combined (e.g. with samtools merge), as are groups 2 and 3. These combined groups are turned into separate BigWig files, representing the "genome-oriented" and "genome-reversed" read pairs.

$\endgroup$

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.