This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Multimapped reads in STAR alignment and subread::featureCounts()

Hello,

My question is where my "Multi-mapped reads" went after STAR Alignment to Quantification with Rsubread featureCounts?

I am aligning my 151 bp untrimmed paired reads to the Human Primary Assembly from Gencode and then quantifying the reads using Rsubread's featureCounts().

I ran STAR with the following parameters:

STAR --genomeDir /storage/genomeDir --runThreadN 6
--readFilesIn /storage/mock1.AAGCATCCTG-GACCGATACA_R1.fastq  /storage/merged/mock1.AAGCATCCTG-GACCGATACA_R2.fastq
--outFileNamePrefix /storage/mock1.AAGCATCCTG-GACCGATACA_
 --outSAMtype BAM SortedByCoordinate --outSAMunmapped Within --outSAMattributes Standard 
--outFilterMatchNminOverLread 0 --outFilterScoreMinOverLread 0 --outFilterMatchNmin  40 --outReadsUnmapped Fastx

I set a cutoff of 40 bp that must align for a read pair to be considered mapped because my reads are untrimmed and the default 66% of the read pair required to map was causing many to be considered "too short" & unmapped. I am not sure of how to appropriately choose the best parameter so any advice of read exploration would be appreciated.

I got 85% uniquely mapped and 13.8% multi-mapped so the results seem good.

I then quantified the pairs using featureCounts in R:

mock2_quant <- featureCounts(files="mock1.AAGCATCCTG-GACCGATACA_Aligned.sortedByCoord.out.bam",isPairedEnd=TRUE, GTF.featureType="gene",GTF.attrType="gene_id",
                       annot.ext= "/storage/Genome_files/gencode.v41.primary_assembly.annotation.gtf" , isGTFAnnotationFile=TRUE)

In my output, I get the following:

Assigned                    20909396
Unassigned_Unmapped         64991
Unassigned_MultiMapping         0
Unassigned_NoFeatures           7240676
Unassigned_Overlapping_Length   0
Unassigned_Ambiguity            6649467

This seems fine but why is unassigned_MultiMapping equal to 0 if 13% from the STAR Alignmnent were seen to be Multimappers? Did I use a parameter that is including these with the other categories in the output such as "Aligned"?

multimap star featurecounts alignment rsubread
• 3,267 views • link
written by Saran • Saran

2 answers

If you do not use argument -M you do not consider multimapping reads in featureCounts. (At least in the command line, never used featureCounts in R tho).

• 0 views • link
written by Marco Pannone • Marco Pannone

I thought I was missing a parameter, thank you!

• 0 views • link
written by Saran • Saran

Hi Marco,

I re-ran featureCounts and added the parameter countMultiMappingReads = TRUE which is equivalent to -M in R but I got the exact same results....any thoughts?

• 0 views • link
written by Saran • Saran

I figured out that if you specify countMultiMappingReads = TRUE with Rsubread::featureCounts() then these counts will be added to your 'Aligned' counts. If you set it to False, then they will display in Unassigned_MultiMapping as shown below after I added the parameter: countMultiMappingReads = FALSE.

Assigned               15124819
Unassigned_Unmapped    64991
Unassigned_MultiMapping    10844204
Unassigned_NoFeatures      4538772
Unassigned_Ambiguity       4291744

Now a very large portion of my reads are seen to be multimappers.... on to the next issue.

• 0 views • link
written by Saran • Saran

UPDATE: The count of fragments mapping to multiple locations is inflated because there is a count for every location that it is mapping; its important not to look at those results as if its a count per read/fragment.

• 0 views • link
written by Saran • Saran

Log in to answer this question.

More posts like this