Hi everybody, I have a set of RNA-seq (single end files). After mapping with STAR (himan genome)
STAR --genomeDir /index --readFilesCommand zcat --readFilesIn file1.fastq.gz --runThreadN 16 --outSAMtype BAM SortedByCoordinate --outWigType bedGraph
I have a huge % of reads mapped to multiple loci
Number of input reads | 49078760
Average input read length | 74
UNIQUE READS:
Uniquely mapped reads number | 9806908
Uniquely mapped reads % | 19.98%
Average mapped length | 73.42
Number of splices: Total | 1022832
Number of splices: Annotated (sjdb) | 985106
Number of splices: GT/AG | 1005417
Number of splices: GC/AG | 6850
Number of splices: AT/AC | 496
Number of splices: Non-canonical | 10069
Mismatch rate per base, % | 0.94%
Deletion rate per base | 0.02%
Deletion average length | 2.33
Insertion rate per base | 0.01%
Insertion average length | 1.22
MULTI-MAPPING READS:
Number of reads mapped to multiple loci | 36919747
% of reads mapped to multiple loci | 75.23%
Number of reads mapped to too many loci | 186914
% of reads mapped to too many loci | 0.38%
UNMAPPED READS:
% of reads unmapped: too many mismatches | 0.00%
% of reads unmapped: too short | 3.97%
% of reads unmapped: other | 0.45%
I would be grateful for any suggestions.
•
written
by
Lila M
1 answer
Reads that map to more than 10 loci are counted as mapping to too many loci. You can change --outFilterMultimapNmax to increase the threshold. Make sure to also increase winAnchorMultimapNmax if you increase it to greater than 50.
I've added the options:
STAR --genomeDir /index --readFilesCommand zcat --readFilesIn file.fastq.gz --runThreadN 16 --outFilterMultimapScoreRange 1 --outFilterMultimapNmax 20 --alignIntronMax 500000 --alignMatesGapMax 1000000 --outSAMtype BAM SortedByCoordinate --outWigType bedGraph
But the results are the same (any other parameter to include?) ... so I will try with BBDuk...what options do I have if those are rRNA?
There are more than 20 (though I don't know how many more) rRNA positions on the assembled human genome. I got rRNA mappings when I set --outFilterMultimapNmax 200, though it is possible you will get them with less than 200. Did you check for rRNA contamination with BBDuk, anyway?
P.S.: I completely agree with genomax.
That said, if you have 75+% rRNA then it may be time to start thinking about doing this experiment over.
yes, I did
bbduk.sh in1=file.fastq.gz outm=ribo.fa outu=nonribo.fa ref=human_ribosomal.fa
been human_ribosomal.fa
and I get:
Input: 49078760 reads 3667028636 bases.
Contaminants: 38486994 reads (78.42%) 2895991336 bases (78.97%)
Total Removed: 38486994 reads (78.42%) 2895991336 bases (78.97%)
Result: 10591766 reads (21.58%) 771037300 bases (21.03%)
So for sure they are rRNA right?
Yes, all of them have a similar contamination level because of the kit (I think). One more question, regarding the command line
bbduk.sh in1=file.fastq.gz outm=ribo.fa outu=nonribo.fa ref=human_ribosomal.fa
which is better, working with the file.fastq.gz (for STAR/salmon) or use directly nonribo.fa to be sure that all the contamination has been removed?
Just to make sure. This isn't an experiment about rRNA correct? Otherwise we would be wrong.
If you have already done STAR alignments you could ignore rRNA counts since the alignment % more or less appears to match bbduk's evaluation. You could use nonribo.fa, if you were going to redo the alignments.
Regarding this answer, apparently STAR has not ignore the rRNA, because after generating the bam file and run samtools view Aligned.sortedByCoord.out.bam | cut -f1 |sort | uniq |wc -l I get a total of 44795631 that I assume are the mapped reads. So STAR is not ignoring at all the rRNA, rigth?
I think so, because I've used the Genome sequence (GRCh38.p10) and Comprehensive gene annotation (gencode.v27.annotation.gtf)
If you use decontaminated files (non_ribo.fq.gz) then you should have only 25% of the original reads in them so the resulting BAM's should be smaller, even when unaligned reads are included. rRNA should not matter now since most reads would not align to it. Am I missing something?
sorry, my mistake. I've run STAR with --outReadsUnmapped in the original fastq.gz file, and the result is the same I get when run STAR without that parameter, so this is why I'm thinking that the genome that I am using for mapping should include rRNA. I hope a much smaller bam file for non_ribo.fq.gz
With samtools flags, filtering options (-f, -F and -G) and the mighty power of internet search engines:
Question: How To Filter Mapped Reads With Samtools
Question: How to extract all the mapped reads from bwa output
Log in to answer this question.
-
Decent Uniquely mapped reads but no features
written by ju_ra •
-
RNA Seq analysis
written by Micro_bioinfo •
-
trouble with 'too short' reads using STAR. each mate maps fine, but the paired end mapping does not
written by james.alan.gregory •
-
STAR mapping results
written by mxlsherry1992 •
-
Terrible RNA mapping result by STAR
written by 610225668 •
-
the lower unique mapping rate after Trimming
written by lilingjoyo •
-
Very low unique mapped reads - STAR alignment
written by jsl •
-
STAR outputs interpretation
written by XBria
-
Long Read Length, yet STAR says many reads too short
written by dec986
-
Problems with STAR mapping output
written by manekineko
First thing to check is if those are rRNA reads.
How can I check it?
Using BBDuk.
You can find entire sequence of human rDNA repeat here. Use it with bbduk.sh.
Was anything suspicious in the FastQC report? Have you checked for rRNA or other contamination?
For RNA-Seq, I'd use the ENCODE settings described in the STAR manual. This parameter set worked well for me in various experiments.