Hi,
I came across the analysis of RNA-seq data, and the trimming of it makes me wonder a bit. We would usually trim the first 10 - 15 base pairs given, the graphic. However, with this, RNA STAR maps only to ~60 % unique reads. If we remove the trimming of the first 10 - 15 base pairs, the unique mapped reads are above 90%! Even if we remove just 1 base pair in the trimming, the quality of the unique mapped reads drops down to 64%. Can someone explain to me why even a removal of 1 basepair changes the result that much? Given how mappers work, with k-mers and suffix arrays, I am surprised the removal of 1 basepair changes the result that much.
Also, why do we have such a non-uniform distribution at the reads for the first 10 - 15 base pairs? What makes them so special that they are actually very important?
Thanks a lot!
3 answers
OK, I think I have it figured out.
Looking at a few reads, the paired end reads are "dovetailing".
Let's say you have 10-bp long reads and you have already trimmed them (by 1 bp) so that they are 9-bp long. The R1 read of TGGTGAGGT and its mate R2 is CCTCACCAG (the reverse complement of R2 is CTGGTGAGG).
That means your alignment would look something like this:
R1: TGGTGAGGT
R2: CTGGTGAGG
Notice the overhang and how the mates extend past each other (i.e. dovetail)? STAR will report this as unaligned.
Now, why does trimming affect these results? Well, in this case, if you hadn't trimmed off that 1 bp, the mates wouldn't extend past each other! (In a lot of cases, I'm seeing these 1 bp overhang results).
I aligned your BAM file of the STAR-unaligned reads with hisat2 and then used samtools to give me the insert size distribution. The vast majority of records had an insert size of 0 (again, this is due to what I described above).
What's the solution? You have to tell STAR to align these by using either the --alignEndsProtrude ConcordantPair option or, perhaps better yet, the --peOverlapNbasesMin option. Using these options, you can specify the allowable overhang or overlap.
In any case, the reason your reads look like this is due to a library prep issue in the lab.
Wow, thank you so much!
For testing, I trimmed now on tail 75 bp away, and the unique read alignment is up to 88%, even if we trim the head 12 too. Which means it confirms your analysis.
What could be a cause for this in the library prep? A too short fragment size in respect to the read length?
What could be a cause for this in the library prep? A too short fragment size in respect to the read length?
If these are FFPE samples then that is normal(ish). It could also be that the starting material was over fragmented and/or the library has not been treated to remove short fragments. Short fragments cluster preferentially.
why do we have such a non-uniform distribution at the reads for the first 10 - 15 base pairs?
Please see the blog post by authors of FastQC located here: https://sequencing.qcfail.com/articles/positional-sequence-bias-in-random-primed-libraries/
If we remove the trimming of the first 10 - 15 base pairs, the unique mapped reads are above 90%!
That is odd. What do you mean by "unique mapped" reads?
I think they mean the uniquely mapped reads that is outputted in the log final at the end of a STAR run.
But that is indeed very odd behavior. I think the best way to diagnose is to literally just use samtools to view the output BAM files and investigate what reads are being mapped in one scenario but unmapped in the other scenario (and see exactly where in the read the mapping is occurring and see if it exists in the genome that was indexed).
Thank you for your reply. I will specify my actions a bit more:
- Trimming with fastp and no hard removal of basepairs:
fastp --thread 8 --report_title 'fastp report for L1_R1_001_fastq_gz.fastq.gz' -i 'L1_R1_001_fastq_gz.fastq.gz' -o first.fastq.gz -I 'L1_R2_001_fastq_gz_R2.fastq.gz' -O second.fastq.gz --detect_adapter_for_pe -f 0 -F 0 - Trimming with fastp and hard removal of 12 basepairs:
fastp --thread 8 --report_title 'fastp report for L1_R1_001_fastq_gz.fastq.gz' -i 'L1_R1_001_fastq_gz.fastq.gz' -o first.fastq.gz -I 'L1_R2_001_fastq_gz_R2.fastq.gz' -O second.fastq.gz --detect_adapter_for_pe -f 12 -F 12
Using RNA STAR for mapping against hg38 with the following commands:
STAR --runThreadN 8 --genomeLoad NoSharedMemory --genomeDir 'hg38_index' --sjdbOverhang 100 --sjdbGTFfile hg38_knownGenes.gtf --sjdbGTFfeatureExon 'exon' --readFilesIn forward.fastq.gz reverse.fastq.gz --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --twopassMode None '' --quantMode - --outSAMattrIHstart 1 --outSAMattributes NH HI AS nM ch --outSAMprimaryFlag OneBestScore --outSAMmapqUnique 60 --outSAMunmapped Within --outBAMsortingThreadN 8 --outBAMsortingBinsN 50 --winAnchorMultimapNmax 50 --limitBAMsortRAM 1000000 && samtools view -b -o mapped.bam Aligned.sortedByCoord.out.bam
This leads in the first case to a unique mapped reads of RNA STAR of Uniquely mapped reads % | 92.26% and in the second case to Uniquely mapped reads % | 54.35%. If we trim just 1 base pair, we get Uniquely mapped reads % | 64.85%.
In the trimming cases, we get e.g. for removal of 12 base pairs we get also a very high reads too short rate Number of reads unmapped: too short | 9856325
% of reads unmapped: too short | 42.06%
For control, we used also TopHat2 on the no hard threshold data vs hard threshold of 12 with the command: tophat2 --num-threads 8 -r 300 --mate-std-dev=20 bowtie2_index/hg38/hg38 forward.fastq.gz reverse.fastq.gz. This leads to 83.9% concordant pair alignment rate. for the non-hard trimming case, and 85.7% concordant pair alignment rate. in case of the 12 base pair removal.
This gives the insight, that trimming the 'noise' seems to work, and improves the mapping for TopHat2. However, it does not give me any idea why RNA STAR acts so differently.
I hope this additional data might give you details to help to explain this behaviour.
Thanks a lot!
Disclaimer: Data was processed with my own Galaxy instance, and I edited the commands for better readability.
Please use hisat2 if you wish to use an alternate aligner for comparison with STAR. Tophat2 is an older aligner that is deprecated.
Note from the TopHat devs says:
Please note that TopHat has entered a low maintenance, low support stage as it is now largely superseded by HISAT2 which provides the same core functionality (i.e. spliced alignment of RNA-Seq reads), in a more accurate and much more efficient way.
As noted by @dsull above you should check the alignments and see why reads are uniquely mapping in one instance and not other.
Thank you for this advice.
Even if Tophat2 is an older aligner and is deprecated, I don't think this gives us an explanation why RNA STAR behaves in the way it does vs. why we don't see this effect with Tophat2.
We want to use RNA STAR in production, and Tophat2 is only referred to as an alternative to investigate if the drop of quality in the alignment comes from the data or the mapper.
RNA STAR parameters: Besides some other output parameters for the SAM/BAM files, the parameters given here are default values of RNA STAR. I confirmed this with the official documentation. https://physiology.med.cornell.edu/faculty/skrabanek/lab/angsd/lecture_notes/STARmanual.pdf
Hisat2:
I used now also hisat2. The following is the result of Trimming with fastp and no hard removal of basepairs as given above:
23460158 reads; of these:
23460158 (100.00%) were paired; of these:
933198 (3.98%) aligned concordantly 0 times
20953931 (89.32%) aligned concordantly exactly 1 time
1573029 (6.71%) aligned concordantly >1 times
----
933198 pairs aligned concordantly 0 times; of these:
37659 (4.04%) aligned discordantly 1 time
----
895539 pairs aligned 0 times concordantly or discordantly; of these:
1791078 mates make up the pairs; of these:
1283064 (71.64%) aligned 0 times
469095 (26.19%) aligned exactly 1 time
38919 (2.17%) aligned >1 times
97.27% overall alignment rate
And this for Trimming with fastp and hard removal of 12 basepairs as given in detail above:
23435619 reads; of these:
23435619 (100.00%) were paired; of these:
11251629 (48.01%) aligned concordantly 0 times
11322054 (48.31%) aligned concordantly exactly 1 time
861936 (3.68%) aligned concordantly >1 times
----
11251629 pairs aligned concordantly 0 times; of these:
7037784 (62.55%) aligned discordantly 1 time
----
4213845 pairs aligned 0 times concordantly or discordantly; of these:
8427690 mates make up the pairs; of these:
1208921 (14.34%) aligned 0 times
5646578 (67.00%) aligned exactly 1 time
1572191 (18.66%) aligned >1 times
97.42% overall alignment rate
The used hisat2 command in version 2.2.1
hisat2 -p 8 -x 'isat2_index/hg38/hg38' -1 'input_f.fastq.gz' -2 'input_r.fastq.gz' --summary-file summary.txt
Hi,
thanks for the offer. I extracted all reads which do no align, and prepared a BAM file with it. https://drive.google.com/drive/folders/181QP93bskLka5WXnKQ1PlvFRykXORNnR?usp=sharing
The data this is based on, is trimmed with fastp with -f 1 and -F 1 and leads to a unique mapping of 65%.
Thanks you for your help.
I expect the issue is result of your unspecified trimming process. You probably didn't really need to trim. It sounds like you are not using options
--clip3pNbases or --clip5pNbases. I think you should consider to do so if you are intent upon dialing in a fixed number of bases to ignore in the alignment.
•
written
by
Malcolm.Cook
Dear Malcolm,
please see my other comment. I added the exact commands I used for the trimming and alignment. In comparison to RNA STAR, TopHat2 does not have any issues with my trimming.
Even though I do not need trimming, I would like to know why even the trimming of one base pair, leads to a drop in alignment rate from 92% to 64%, trimming 12 base pairs leading to only 54% alignment rate. At the same time TopHat has 83% for non-trimmed and 85% for 12 base pair trimmed data.
Why is RNA STAR behaving like this?
Thanks a lot!
Log in to answer this question.
-
Merge paired end reads after mapping in BAM files
written by muhammad_elhossary •
-
Read length for single -ended RNA seq data analysis
written by vinaya •
-
STAR aligner can't map too short reads
written by Assa Yeroslaviz
-
Extract mapped SE and properly paired PE reads from a bam file
written by Jautis
-
ATAC-seq problem with STAR - Multimapping reads
written by Rafael Soler •
-
filter internal priming events in oligdT based Reverse transcription of 3'end library
written by linjc.xmu •
-
the lower unique mapping rate after Trimming
written by lilingjoyo •
-
How to extract uniquely mapped reads in Bowtie2
written by Ankit
-
Different trimming SLIDINGWINDOW threshold value mapping comparison
written by Hughie •
-
Recommended doc for interpreting STAR alignment stats
written by emblake •
Could be helpful if you provide the STAR output summary (e.g. after trimming, what do these reads then get classified as?)
See my other comment, but also the full statistics: STAR output for removingno base pair with fastp (-f 0 and -F 0):
STAR output for removing one 1 base pair with fastp (-f 1 and -F 1):
STAR output for removing 12 base pair with fastp (-f 12 and -F 12):
Output for TopHat2 on non-removal:
Output for TopHat2 on 12 base-pair removal: