HTSeq analysis for Illumina tissue samples?

Question

[pbio]

Hi,

I have downloaded a few tissue bam files from illumina body map from the link ftp://ftp.ensembl.org/pub/release-81/bamcov/homo_sapiens/genebuild/

Now to get the read counts from each of the tissue sample bam files, I did HTSeq analysis. But unfortunately for a few tissues it worked really well but for a few tissues it ended up showing different error each time I run HTSeq analysis.

The HTSeq command I am using for processing bam files is

samtools view GRCh38.illumina.ovary.1.bam |\
  htseq-count -i gene_id - Homo_sapiens.GRCh38.80.gtf > htseq.ovary.results.txt

shows an error

Error occured when processing SAM input (line 2):
  'pair_alignments' needs a sequence of paired-end alignments
  [Exception type: ValueError, raised in __init__.py:603]

I tried sorting the sam files too, but still it didn't work :( Showed error:

Error occured when reading beginning of SAM/BAM file.
  [Exception type: StopIteration, raised in count.py:88]

How to solve this issue please help, someone!

Answer 1

[GouthamAtla]

Error at the end of HT-seq Analysis and output.txt with 0 bytes in size

[pbio]

Error occured when reading beginning of SAM/BAM file.
  ('SAM line does not contain at least 11 tab-delimited fields.', 'line 1')
  [Exception type: ValueError, raised in _HTSeq.pyx:1276]

[GouthamAtla]

The error is clear this time. SAM line does not contain at least 11 tab-delimited fields.', 'line 1' there is something wrong with your file. Why don't you paste few lines

samtools view <in.bam> | head

[pbio]

I am sorry for the delay in response, this is how to file looks like

HWI-BRUNOP16X_0001:3:64:13105:149752#0  97      KI270757.1      49010   0       50M     8       89223567        0       GCCAAAATTGACAAAT               GGGATCTAATTAAACTAAAGAGCTTCTGCACAGT      eaceaddfffdggggdggcdgfgfgggggfggggfdgggdgggdgggggh      RG:Z:ovary_50_fca       XT:A:R  NM:i:0 S               M:i:0   AM:i:0  X0:i:69 XM:i:0  XO:i:0  XG:i:0  MD:Z:50
HWI-BRUNOP16X_0001:6:25:9464:120176#0   0       KI270757.1      63302   25      75M     *       0       0       TGGTGAGCAATTTTTTCATGTGTT               TTTTGATTTGCATTTCTCTGATGGCCAGTGATGGTGAGCAATTTTTTCATG     gaggggggggggggggggggggggggfffgggggggggggfgggeggggggggggggghgggggegggggggggg    R               G:Z:ovary_75_fcb        XT:A:U  NM:i:37 X0:i:1  X1:i:0  XM:i:37 XO:i:0  XG:i:0  MD:Z:9T20C1G0C0A0T0G0A0A0T0G0T0C0T0T0C0T0T0T0T0G0A0G0A1G               0T2C0T1T0T0C0A1G1C0C0T0T2
HWI-BRUNOP16X_0001:6:66:18737:92611#0   0       KI270757.1      67475   0       75M     *       0       0       GATGTGGAGAAATAGGAACACTTT               TACACTGTTGGTGGGACTGTAAACTAGTTCAACCATTGTGGAAGTCAGTGT     ggggggggghgggggggggggggggfgggggggggfgggggggfgggggggggghgggggggefggggggggfef    R               G:Z:ovary_75_fcb        XT:A:R  NM:i:0  X0:i:4670       XM:i:0  XO:i:0  XG:i:0  MD:Z:75
HWI-BRUNOP16X_0001:6:43:20784:119867#0  0       KI270757.1      67490   0       75M     *       0       0       GAACACTTTTACAGTGTTGGTGGG               ACTGTAAACTAGTTCAACCATTGTGGAAGTCAGTGTGGCGATTCCTCAGTG     gfggggggggggggggggggegggggggggggggggggggfgggfggggggggggfgfgfegggaeggggggghe    R         

[pbio]

Thanks, the issue was solved using the following commands

samtools view -bf 1 GRCh38.illumina.ovary.1.bam > ovary.bam
samtools view ovary.bam | python -m HTSeq.scripts.count --type=CDS --idattr=gene_id --mode=union --stranded=yes - Homo_sapiens.GRCh38.80.gtf