Hello! I am analyzing smallRNA data of human genome. I have aligned my read against genome with bowtie2:
bowtie2 --end-to-end --very-sensitive -x hg38 -U input_file -S output_file_align
and I used featurescount in order to obtained the counts. For the counting: I used the gtf file with miRNA: featureCounts -T 5 -t exon -g gene_name -a GTF_miRNA -R SAM -o fileOutput output_file_align and the reads not assigned to miRNA were annotated with the gtf from gencode for the other non-coding RNA molecules.
This is the result of the alignment:
41798971 reads; of these:
41798971 (100.00%) were unpaired; of these:
1101846 (2.64%) aligned 0 times
10724437 (25.66%) aligned exactly 1 time
29972688 (71.71%) aligned >1 times
97.36% overall alignment rate
This is the results of featurescount for miRNA:
Total reads : 41798971
Successfully assigned reads : 15157113 (36.3%)
And this for the other non-coding RNA molecules:
Total reads : 41798971
Successfully assigned reads : 12408845 (29.7%)
As reported in the featurescount manual:
For the reporting of a multi-mapping read, users may choose to not report any alignments for the read (by default) or report up to a pre-defined number of alignments (‘–multiMapping’ and ‘-B’ options).
How is it possible that I obtained the previously result being the 71.71% of read aligned >1 times?
Thank you in advance
1 answer
reads not assigned to miRNA were annotated
I don't think that's right since the number of total reads is the same in both featureCount cases
In regard to featureCounts not recognizing the multi-mappers -- I'm not sure how what kind of information featureCounts is relying on (re-check the manual, maybe it mentions it). It may need the flag that is typically set by Picards MarkDuplicates, which bowtie2 probably does not do (if I recall correctly, bowtie2 reports multi-mappers via optional tags in the SAM file (e.g. XS), not via the flag in the fifth column of the SAM file).
Log in to answer this question.