This is a test version of Biostars. For the public version, visit https://www.biostars.org.
HTseq error with paired end alignments

I am running HTseq code as follows:

python2.7 -m HTSeq.scripts.count -f bam -r pos -s no /extstor/sudhirlab/Reety/Bhavana/Set3_firstAgrigenome/Sorted_CFPAC-1_N1_aligned.bam /extstor/sudhirlab/Reety/Bhavana/Indexed_Genomes/GRCh38/ChromosomeNamed_Gtf_GRCh38.gtf > /extstor/sudhirlab/Reety/Bhavana/Set3_firstAgrigenome/Counts_Sorted_CFPAC-1_N1_aligned.txt

When I execute the command, the command starts to run showing the following:

100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
500000 GFF lines processed.

....and so on

However, at the end it is giving me the following warning:

Warning: Mate records missing for 153562 records; first such record: <SAM_Alignment object: Paired-end read 'HWI-1KL120:151:C3350ACXX:6:2205:21118:19311' aligned to chrMT:[5593,5692)/+>.
34600000 SAM alignment record pairs processed.
Warning: Mate pairing was ambiguous for 133693 records; mate key for first such record: ('HWI-1KL120:151:C3350ACXX:6:2301:10892:62181', 'second', 'chr1', 14614, 'chr1', 14800, 285).
34697807 SAM alignment pairs processed.

When I open my output text file, I see that there are several rows of:

ENSG00000000003 0
ENSG00000000005 0
ENSG00000000419 1504
ENSG00000000457 409
ENSG00000000460 224

And in the end this is what it says:

__no_feature    3952696
__ambiguous     2186111
__too_low_aQual 1112784
__not_aligned   235455
__alignment_not_unique  1318607

So has my HTseq worked or not? I am very confused. I am guessing these warning signs mean that HTseq cannot handle any unpaired reads? How do I proceed with this?

Second, when I am running the same on another set of sorted bam files, then I am getting a different error. This second set of sorted bam files, originally were created from sam files generated after alignment with Hisat2, wherein I had input both paired and unpaired outputs (given by trimmomatic), to use for performing alignment. The error is as follows:

Error occured when processing SAM input (record #10 in file /extstor/sudhirlab/Reety/Bhavana/Set3_firstAgrigenome/Sorted_OZ-GCTM5-POSITIVE_aligned.bam):
  Sequence of paired-end alignments expected, but got single-end alignment.
  [Exception type: ValueError, raised in __init__.py:761]

And the corresponding .txt file in which I had directed my output to, is empty.

My guess after reading some other posts is that HTseq does not handle unpaired reads well. So what can I do about it ? I am very new to informatics, so any help will be great.

rna seq rna-seq htseq
• 393 views • link
written by bnayer26 • bnayer26

I think HTSeq has no problem with single-end reads per se, the problem is mixing single-end and paired-end in the same file. Even if HTSeq didn't have a problem with this, you shouldn't mix single- and paired-end counts, because you may have a batch effect, due to the different "mappabilities" of each type of read.

• 0 views • link
written by h.mon h.mon

I don't know at what step this "mixing" of single end reads would have happened? Could you elaborate. My raw data consisted of paired end reads only. I used that for adapter trimming with Trimmomatic, which gave me a "paired_trimmed" and an "unpaired trimmed" file for each read. Then I used both paired trimmed and unpaired trimmed for Read 1 and Read 2 (total 4 files) and used them in the alignment command with Hisat2. The output was a sam file for each sample, which i converted to bam with samtools view command then i used samtools sort to sort by coordinate (default) and then used this file for htseq. So i am not sure what to change in this pipeline

• 0 views • link
written by bnayer26 • bnayer26

0 answers

No answers yet.

Log in to answer this question.

More posts like this