HTSeq-Count: no_feature too high?

Question

[sea.joson]

Hello!

I am learning how to do differential expression. I was given transcriptome data and I am fortunate that the sample already has an established annotated genome. So far, I aligned the transcripts to the reference genome via STAR and I've had good mapping stats (>95%). Now, I'm quantifying the counts with HTSeq-Count. I expected it to be a relatively smooth process, but the results have been bugging me.

Here is the line I used:

htseq-count -f bam -r pos -t exon $RNA/AN0_1.Aligned.sortedByCoord.out.bam $GEN > AN0_1.no.exon.counts

The results:

__no_feature    14805713
__ambiguous     4766268
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  113292

I feel like that's a high value for no_features right?

I've tried playing around with the strandedness (-s), but all options give similar results. I also went back and made sure the GTF file has the gene_id attribute. I also tried using a GFF file instead.

Might it also have to do with this warning sign from the err file?

Warning: Mate pairing was ambiguous for 23164 records; mate key for first such record: ('A00564:87:HJKFCDSXX:2:2205:13838:9330', 'first', 'CP059858.1', 91661, 'CP059858.1', 91831, 317).

Any help will be much appreciated!

[Marco Pannone]

Warning: Mate pairing was ambiguous

It makes me think you have paired-end reads that have been trimmed inappropriately without using paired-end mode.

Also, can you try running it again with -t gene and let me know what results you get?

[sea.joson]

This is run with -t gene

__no_feature    14786206
__ambiguous     4691
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  113292

*Warning: Mate pairing was ambiguous for 23164 records; mate key for first such record: ('A00564:87:HJKFCDSXX:2:2205:13838:9330', 'first', 'CP059858.1', 91661, 'CP059858.1', 91831, 317).*

I'm not super certain about the trimming because when I got the files, I think they've already been pre-processed a bit. I checked them out using fastqc and they seem fine.

https://imgur.com/n6Zyg1q https://imgur.com/h6JcYzP

[GenoMax]

Do you have identical number of reads in R1/R2 files?

Can you run the following program from BBMap suite on each file and show the output

$ readlength.sh -Xmx2g in=your_file.fastq.gz
Executing jgi.MakeLengthHistogram [-Xmx2g, in=your_file.fastq.gz]


Processed 36129 reads.
#Reads: 36129
#Bases: 10838700
#Max:   300
#Min:   300
#Avg:   300.0
#Median:        300
#Mode:  300
#Std_Dev:       0.0

[sea.joson]

readlength.sh in=AN0_1_1.fq.gz in2=AN0_1_2.fq.gz

Processed 56800750 reads.
#Reads: 56800750
#Bases: 8520112500
#Max:   150
#Min:   150
#Avg:   150.0
#Median:        150
#Mode:  150
#Std_Dev:       0.0
#Read Length Histogram:
#Length reads   pct_reads       cum_reads       cum_pct_reads   bases   pct_bases       cum_bases       cum_pct_bases
150     56800750        100.000%        56800750        100.000%        8520112500      100.000%        8520112500      100.000%

[GenoMax]

Can you run the command independently on the two files?

[sea.joson]

readlength.sh in=AN0_1_1.fq.gz

Processed 28400375 reads.
#Reads: 28400375
#Bases: 4260056250
#Max:   150
#Min:   150
#Avg:   150.0
#Median:        150
#Mode:  150
#Std_Dev:       0.0
#Read Length Histogram:
#Length reads   pct_reads       cum_reads       cum_pct_reads   bases   pct_bases       cum_bases       cum_pct_bases
150     28400375        100.000%        28400375        100.000%        4260056250      100.000%        4260056250      100.000%
Time:   67.509 seconds.

readlength.sh in=AN0_1_2.fq.gz

Processed 28400375 reads.
#Reads: 28400375
#Bases: 4260056250
#Max:   150
#Min:   150
#Avg:   150.0
#Median:        150
#Mode:  150
#Std_Dev:       0.0
#Read Length Histogram:
#Length reads   pct_reads       cum_reads       cum_pct_reads   bases   pct_bases       cum_bases       cum_pct_bases
150     28400375        100.000%        28400375        100.000%        4260056250      100.000%        4260056250      100.000%
Time:   69.977 seconds.

[GenoMax]

These files appear to be non-trimmed and intact.

Did you make sure that these actually have reads from R1 and R2 by inspecting the headers? Hopefully you don't have the same data in the two files.

Also inspect the aligned BAM in IGV to see where the reads are aligning (i.e. are they piling up on exons etc). Do you have rRNA gene models in your GTF? Perhaps you may have a lot of rRNA contamination.

[sea.joson]

Oh. I did not trim my reads. Do you recommend I try again but with trimmed reads instead? Also, may I ask how you are able to tell that the files "appear to be non-trimmed and intact?"

I did check the R1 and R2 and they are not the same data.

I also went ahead and checked the bam file in IGV. Unfortunately, I do no yet know how to fully interpret what I was seeing.

[GenoMax]

Do you recommend I try again but with trimmed reads instead?

No need. STAR should be able to handle any adapter sequences present by soft-clipping them.

Also, may I ask how you are able to tell that the files "appear to be non-trimmed and intact?"

You can see that the length of reads is identical (150 bp ) in both files and number of bases is also identical. If reads are trimmed to remove adapters there will be a range of sizes present and the number of bases will likely be different in two files.

I do no yet know how to fully interpret what I was seeing.

Find a smallish gene in IGV (which shows exon structure) and post a screen shot of what the alignments look like.

[sea.joson]

Is this right? https://imgur.com/hXGSf4b

[GenoMax]

That looks right. Reads piling up under exons and a good number too.

If you are open to trying another counting program out can you try featureCounts? It is part of a package called subread (LINK). A tutorial is available: https://subread.sourceforge.net/featureCounts.html

[sea.joson]

Yeah I think this is a next move for me.

[sea.joson]

Finally got to try featureCounts. It seems to be better, correct?

Status  AN0_1.Aligned.sortedByCoord.out.bam
Assigned        45141493
Unassigned_Unmapped     0
Unassigned_Read_Type    0
Unassigned_Singleton    0
Unassigned_MappingQuality       0
Unassigned_Chimera      0
Unassigned_FragmentLength       0
Unassigned_Duplicate    0
Unassigned_MultiMapping 226584
Unassigned_Secondary    0
Unassigned_NonSplit     0
Unassigned_NoFeatures   8469135
Unassigned_Overlapping_Length   0
Unassigned_Ambiguity    24850

[sea.joson]

https://imgur.com/fxxd4L2

[Marco Pannone]

40M reads assigned is great.

[sea.joson]

Thank you both for all the help! Now I just gotta figure why that was the case in the first place.

[Marco Pannone]

From where did you download the .gtf file?

[GenoMax]

sea.joson: And to add another point to this question, do the chromosome names in your reference and GTF file match perfectly.

[sea.joson]

Yep. I checked and made sure that the chromosome names used are the same across the board.

[sea.joson]

https://www.ncbi.nlm.nih.gov/genome/360?genome_assembly_id=968151

Fortunately, the genome has already been worked on. The samples I'm working on are transcriptomes of the same strain grown in different conditions.

[Marco Pannone]

Would you give a try with another .gtf for the same organism from another source? Maybe Ensembl?