Very low successfully assigned alignments with feature counts
Hello everyone,
I am stuck trying to analyze some single-end RNAseq data from human tissue. My issue is that the alignment with HISAT 2 went very well: 94.95% overall alignment rate.
However, when I use featureCounts, I get:
- 5.7% when I set the strandSpecific parameter to 1.
- 5.3% when I set the strandSpecific parameter to 2 (I guess this would be the correct one)
- 18.8% when I set the strandSpecific parameter to 0
Since this is the first time I'm analyzing sequencing data, I was wondering if someone could guide me on what might be happening or if this is normal.
Thank you very much in advance.
I include more information about the kits and commands used:
The kits that have been used:
- For library preparation: Illumina® Stranded Total RNA Prep, Ligation with Ribo-Zero Plus (16 Samples) (REF: 20040525)
- For sequencing: NovaSeq 6000 SP Reagent Kit v1.5 (100 cycles) (REF: 20028401)
In our case, the sequencing was single (R1).
For the alignment, I use HISAT2 with the following command:
hisat2 -q \
--rna-strandness F \
-x /home/manu/research/genome/grch38/genome \
-U /home/manu/research/analisis/data/ITK2-004-S1.fastq.gz | \
samtools sort -o /home/manu/research/analisis/align/demo_trimmed.bam
99524981 reads; of these:
99524981 (100.00%) were unpaired; of these:
5030213 (5.05%) aligned 0 times
58319656 (58.60%) aligned exactly 1 time
36175112 (36.35%) aligned >1 times
94.95% overall alignment rate
Later I use featureCounts to count the reads:
counts <- featureCounts(files = "demo_trimmed.bam",
annot.ext = "Homo_sapiens.GRCh38.110.gtf",
isPairedEnd = FALSE,
strandSpecific = 1,
isGTFAnnotationFile = TRUE,
GTF.featureType = "exon",
GTF.attrType = "gene_id")
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
Rsubread 2.14.2
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| ||
|| demo_trimmed.bam ||
|| ||
|| Paired-end : no ||
|| Count read pairs : no ||
|| Annotation : Homo_sapiens.GRCh38.110.gtf (GTF) ||
|| Dir for temp files : . ||
|| Threads : 1 ||
|| Level : meta-feature level ||
|| Multimapping reads : counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
\\============================================================================//
//================================= Running ==================================\\
|| ||
|| Load annotation file Homo_sapiens.GRCh38.110.gtf ... ||
|| Features : 1649677 ||
|| Meta-features : 62754 ||
|| Chromosomes/contigs : 47 ||
|| ||
|| Process BAM file demo_trimmed.bam... ||
|| Strand specific : stranded ||
|| Single-end reads are included. ||
|| Total alignments : 225796353 ||
|| Successfully assigned alignments : 12844679 (5.7%) ||
|| Running time : 3.32 minutes ||
|| ||
|| Write the final count table. ||
|| Write the read assignment summary. ||
|| ||
\\============================================================================//
\\============================================================================//
counts <- featureCounts(files = "demo_trimmed.bam",
annot.ext = "Homo_sapiens.GRCh38.110.gtf",
isPairedEnd = FALSE,
strandSpecific = 2,
isGTFAnnotationFile = TRUE)
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
Rsubread 2.14.2
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| ||
|| demo_trimmed.bam ||
|| ||
|| Paired-end : no ||
|| Count read pairs : no ||
|| Annotation : Homo_sapiens.GRCh38.110.gtf (GTF) ||
|| Dir for temp files : . ||
|| Threads : 1 ||
|| Level : meta-feature level ||
|| Multimapping reads : counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
\\============================================================================//
//================================= Running ==================================\\
|| ||
|| Load annotation file Homo_sapiens.GRCh38.110.gtf ... ||
|| Features : 1649677 ||
|| Meta-features : 62754 ||
|| Chromosomes/contigs : 47 ||
|| ||
|| Process BAM file demo_trimmed.bam... ||
|| Strand specific : reversely stranded ||
|| Single-end reads are included. ||
|| Total alignments : 225796353 ||
|| Successfully assigned alignments : 34526662 (15.3%) ||
|| Running time : 3.47 minutes ||
|| ||
|| Write the final count table. ||
|| Write the read assignment summary. ||
|| ||
\\============================================================================//
counts <- featureCounts(files = "demo_trimmed.bam",
annot.ext = "Homo_sapiens.GRCh38.110.gtf",
isPairedEnd = FALSE,
strandSpecific = 0,
isGTFAnnotationFile = TRUE)
========== _____ _ _ ____ _____ ______ _____
===== / ____| | | | _ \| __ \| ____| /\ | __ \
===== | (___ | | | | |_) | |__) | |__ / \ | | | |
==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
Rsubread 2.14.2
//========================== featureCounts setting ===========================\\
|| ||
|| Input files : 1 BAM file ||
|| ||
|| demo_trimmed.bam ||
|| ||
|| Paired-end : no ||
|| Count read pairs : no ||
|| Annotation : Homo_sapiens.GRCh38.110.gtf (GTF) ||
|| Dir for temp files : . ||
|| Threads : 1 ||
|| Level : meta-feature level ||
|| Multimapping reads : counted ||
|| Multi-overlapping reads : not counted ||
|| Min overlapping bases : 1 ||
|| ||
\\============================================================================//
//================================= Running ==================================\\
|| ||
|| Load annotation file Homo_sapiens.GRCh38.110.gtf ... ||
|| Features : 1649677 ||
|| Meta-features : 62754 ||
|| Chromosomes/contigs : 47 ||
|| ||
|| Process BAM file demo_trimmed.bam... ||
|| Single-end reads are included. ||
|| Total alignments : 225796353 ||
|| Successfully assigned alignments : 42545035 (18.8%) ||
|| Running time : 3.36 minutes ||
|| ||
|| Write the final count table. ||
|| Write the read assignment summary. ||
|| ||
[\\============================================================================//]\\============================================================================//
• 1,523 views
•
link
updated
by
Istvan Albert
•
written
by
Manuel •
•
written
by
Manuel •
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Mouse genome featureCounts
written by rajdeepboral00 •
-
High Amount of UnAssignedNoFeatures in FeatureCounts Output
written by Tom Brady •
-
featureCounts results: low rate of 'Successfully assigned alignments'
written by Hamtaro •
-
featureCounts Output "The reads are assigned on the single-end mode"
written by ExtentHonest56 •
-
Rsubread featurecounts
written by elengss •
-
Rsubread FeatureCounts return 0.0% assigned
written by Bonita •
-
FeatureCounts dont read all features
written by David •
-
R featureCounts very low assigned read rates
written by 1111111 •
-
Can't find output after featureCounts
written by ovariohisterectomia •
-
FeatureCounts on HG38
written by elb
Is the genome/annotation from the same source i.e. Ensembl or so it appears?
Load up your data in IGV and visualize it relative to your annotation file. It will tell you where the data aligns.
5% match on a 99% alignment is not the expected outcome
it could be that the annotations don't match the genome - that would be my first guess; make sure the chromosome names match, too
(that being said IGV will convert the chr1 to 1 but featureCounts won't)