Junction counting and filter out low mapping quality reads

[Мария Романова]

Hello. I am interested in using featureCounts as a tool for detecting abnormal splicing events in my RNA-seq data. I have a sample with a known mutation which lead to exon skipping and I would have liked to count how many reads span this junction. I am trying now to decrease the number of artifacts i.e junctions that are supported by reads with low mapping quality. For that purpose, I add -Q flag to the command, but I get not what I expected.

My command is:

featureCounts -J -Q 255 -G GRCh38.primary_assembly.genome.fa -a gencode.v32.primary_assembly.annotation.gtf sample.sam -o sample_fc

A one line example from output sample_fc.jcount:

PrimaryGene     SecondaryGenes  Site1_chr       Site1_location  Site1_strand    Site2_chr       Site2_location  Site2_strand    sample.sam

ENSG00000108821.13      NA      chr17   50185429        +       chr17   50186473        +       3308
ENSG00000108821.13      NA      chr17   50185429        NA      chr17   50194770        NA      1
ENSG00000108821.13      NA      chr17   50185429        +       chr17   50298247        +       1
ENSG00000108821.13      NA      chr17   50185429        +       chr17   50399843        +       1
ENSG00000108821.13      NA      chr17   50185429        +       chr17   50425454        +       1

The same command whithout -Q flag outputs the same result.

But when I use samtools view before this with -q 255 flag I got this:

PrimaryGene     SecondaryGenes  Site1_chr       Site1_location  Site1_strand    Site2_chr       Site2_location  Site2_strand    sample.sam

ENSG00000108821.13      NA      chr17   50185429        +       chr17   50186473        +       4

It is obvious that reads are not filtered out for junction counting with featureCounts, thus I wonder if it is not supposed to or I did something incorrectly.

Thank you in advance for your response. A subsample of my sample is in the attachment in case needed

[Yang LIAO]

When counting junctions, featureCounts doesn't apply filters on reads. It uses all mapped reads in the input SAM or BAM files as long as they contain junction(s). 

You may need to filter the reads using samtools and run featureCounts on the results.

Cheers,

Yang

[Wei Shi]

We will also try to support read filtering for junction counting in featureCounts.

[Мария Романова]

I understand. Thank you very much for your response.

пятница, 30 октября 2020 г. в 01:42:35 UTC+3, shiw...@gmail.com: